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FOREWORD 


This report presents results of the expansion and improvement of the 
FORMA system for response and load analysis. The acronym FORMA stands 
for FORTRAN Matrix Analysis. The study, performed from 16 May 1975 
through 17 May 1976 was conducted by the Analytical Mechanics Department, 
Martin Marietta Corporation, Denver Division, under the contract NAS8- 
3137 j. The program was administered by the National Aeronautics and 
Space Administration, George C. Marshall Space Flight Center, Huntsville, 
Alabama under the direction of Dr. John R. Admire, Structural Dynamics 
Division, Systems Dynamics Laboratory* 

This report is published in seven volumes: 

Volume I - Programming Manual, 

Volume IIA - Listings, Dense FORMA Subroutines, 

Volume IIB - Listings, Sparse FORMA Subroutines, 

Volume IIC - Listings, Finite Element FORMA Subroutines, 

Volume ILIA - Explanations, Dense FORMA Subroutines, 

Volume IIIB - Explanations, Sparse FORMA Subroutines, and 
Volume IIIC - Explanations, Finite Element FORMA Subroutines- 
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ABSTRACT 


This report presents techniques for the solution of structural 
dynamic systems on an electronic digital computer using FORMA (FORT RAN 
Matrix Analysis). 

FORMA is a library of subroutines coded in FORTRAN IV for the effi- 
cient solution of structural dynamics problems. These subroutines are 
in the form of building blocks that can be put together to solve a large 
variety of structural dynamics problems. The obvious advantage of the 
building block approach is that pr gramming and checkout time are limi- 
ted to that required for putting the blocks together in the proper order 

The FORMA method has advantageous features such as: 

1. subroutines in the library have been used extensively for many 
years and as a result are well checked out and debugged; 

2. method will work on any computer with a FORTRAN IV compiler; 

3. incorporation of new subroutines is no problem; 

4. basic FORTRAN statements may be used to give extreme flexi- 
bility in writing a program. 

Two programming techniques are used in FORMA: dense and sparse. 
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I. INTRODUCTION 


A listing of the source deck of each finite element FORMA subroutine 
is given in this volume to remove the "black box" aura of the subroutines 
so that the analyst may better understand the detailed operations of each 
subroutine. 

The FORTRAN IV progranming language is used in all finite element 
FORMA subroutines. 
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II. SUBROUTINE LISTINGS 


The subroutines are given in alphabetical order with numbers 
coming before letters. 
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AXIAL — 1/ 4 


SUBROUTINE AXIAL ( XY 2, JDOF ,EUL ,NUTFL ,N J , 

* NUTMX, NUTKX.NUTLT, NUTST, 

* W*T,$,KX,KJ,KE,KW) 

DIMENSION XYZ (KX , 1 ) , JDOFI KJ ,1) *EUL ( KF. ,1 ) ,W( KW , 1) , TIKW ,1 > ,S (KW ,1 ) 
DIMENSION CJ(3,2), EJ(3,21, IV1(6) 

DATA NAME L/6HAX I AL /, NRW ,NRLT/6, 2/ , IBLNK/6H /, KCJ/3/ 

DATA NIT,NCT/5»6/ 

SUBROUTINE TC CALCULATE (UN OPTION) FINITE ELEMENT ... 

MASS MATRICES AND IVECS (ON NUTMX ) , 

STIFFNESS MATRICES (SAME AS GLOBAL LOAD TRANSFORMATION MATRICES) 

AND IVECS (ON NUTKX)» 

LOCAL LOAD TRANSFORMATION MATFICFS AND IVFCS (ON NUTLT), 

STRESS TRANSFORMATION MATRICES AND IVECS (ON NUTST) 

FOR AXIAL *OD ELEMENTS. 

MASS, ST T FFNFSS MATRICFS ARE IN GLOBAL COORDINATE DIRECTIONS. 

GLOBAL COORDINATE ORDER FOR EACH ELEMENT IS 
(U,V,W) JOINT It THEN JOINT 2. 

WHERE U » V , W APE TRANSLATIONS. 

IVEC GIVES ELEMENT DOF INTO GLOBAL DOF. EXAMPLES... 

I VE C ( 6 ) -834- PLACES ELEMENT DOF 6 INTO GLOBAL DOF 834. 

IVEC ( 3 ) =0 OMITS ELEMENT DOF 3 FROM GLOBAL DOF. THIS CONSTRAINS 

ELEMENT DOF 3 TO ZERO MOTION. 

GLOBAL LOAD TRANSFORMATION MATRIX RFLATFS LOADS AT ROD ENDS IN GLOBAL 
COORDINATE DIRECTIONS TO DEFLECTIONS IN THE GLOBAL COORDINATE 
DIRECTIONS. 

ROW ORDER IN GLOBAL LOAD TRANSFORMATION MATRIX IS 
(PU,PV,PW) JOINT 1, THEN JOINT 2. 

WHERE P IS FORCE. 

LOCAL LOAD TRANSFORMATION MATRIX RELATES LOADS AT PCD ENDS IN LOCAL 
COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE 
DIRECTIONS. 

ROW ORDER IN LOCAL LOAD TRANSFORMATION MATRIX IS 
PX1,PX? 

WHEPE PX IS AXIAL FORCE. 

PXl(-), PX2(+) IS TENSION. PX1(+), PX2 (-) IS COMPRESSION. 

STRESS TRANSFORMATION MATRIX RELATES STRESS AT ROD ENDS IN LOCAL 
COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE DIRECTIONS. 

ROW ORDER IN STRESS TRANSFORMATION MATRIX IS 
SIGMA— XI , S1GMA-X2 
WHERE SIGMA IS NORMAL STRESS. 

SXl(-), SX2(+> IS TENSION. SX1UJ, SX2 (-) IS COMPRESSION. 

DATA ARRANGEMENT ON NUTMX, NUTKX » NUTLT, NUTST FOR EACH' FINITE 
ELEMENT IS ( W-M»K ,LT ,ST) 

WRITE (NUTWX) NAMEW ,N EL ,NR, NC ,NAMEL , ( IB LNK, 1 = 1 , 5 ) , 

( iW(I,J),I=l,NR),J=l,NC),(IVEC(I),I=l,NC) 

CALLS FORMA SUBROUTINES MAS1A, PAGFHD, STF1A, ZZBOMB . 

DEVELOPED BY WA EENFIELD, OS BODLEY, RL WOHLEN. JANUARY 1973. 

LAST REVISION BY WA EENFIELD. MARCH 1976. 

****************************************** ************************ ******** 

INPUT DATA READ IN THIS SUBROUTINE FROM NUTEL. IF NUTEL = 5, DATA IS 
READ FROM CAPDS. 

NAMEM, NAM EK,NAMELT, NAME ST 
RC,E 


FORMAT (4 ( A6,4X ) 
FOPMAT <2(5X,E10n 
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20 NEL, Jit J2 ,A1,A2 FORMAT (3I5,2E10) 

IF (J1 .FO. 0) RETURN 
GO TC 20 

DEFINITION OF INPUT VARIABLES. 

NAMEM = TYPF OF MASS MATRIX WANTED. 

= Ml, DIAGONAL LUMPED. 

= M2, CONSISTENT. 

- EH OR 6HNUMASS, NO MASS MATRIX CALCULATED. 

NAMEK = TYPF CF STIFFNESS MATRIX WANTED. 

= Kl, CONSTANT AXIAL FORCE ASSUMED. 

= 6H OR 6HN0STIF, NO STIFFNESS MATRIX CALCULATED. 

MAMELT = IDENTIFICATION NAME FOR LOAD TRANSFORMATION MATRICES. 

= 6H OR 6HN0L0AD, NO LOAD TRANSFORMATIONS CALCULATED. 

NAMEST = IDENTIFICATION NAME FOR STRESS TRANSFORMATION MATRICES. 

= 6H OR 6HN0STRS, NO STRESS TRANSFORMATIONS CALCULATED. 

RC = MASS DENSITY. 

E = YOUNGS MODULUS OF ELASTICITY. 

NEL = FINITE ELEMENT NUMBER. FOP REFERENCE ONLY, NOT USED IN 
CALCULATIONS. WRITTEN ON NUTMX, ETC . 

J1 = JOINT NUMBER AT POD END 1. 

J2 = JOINT NUMBER AT ROD END 2. 

Al = CROSS-SECTION AREA AT ROD END 1. 

A2 - CROSS-SECTION AREA AT ROD END : . 

EXPLANATION OF INPUT FORMATS. NUMBER INDICATES CARD COLUMNS USED. 

I = INTEGER DATA, RIGHT ADJUSTED. 

E = DECIMAL POINT DATA, ANYWHERE IN FIELD. EXPONENT RIGHT ADJUSTED. 
X = CARD COLUMNS SKIPPED. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

XY2 = MATRIX OF JOINT GLOBAL X,Y,Z LOCATIONS. ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 
X , Y , 2 LOCATIONS RESPECTIVELY. SIZE(NJ,3). 

JDOF = MATRIX OF JOINT GLOBAL DEGREES OF FREEDOM. ROWS CORRESPOND 
TC JOINT NUMBFRS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 
TRANSLATION DOFS AND COLUMNS 4,5,6 CORRESPOND TO THE JOINT 
ROTATION DOFS. SIZE(NJ,6>. 

EUL = MATRIX OF JOINT EULER ANGLES (DEGREES). ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE 
GLOBAL X , Y, 2 PERMUTATION. S1ZF(NJ,3). 

NUTEL = LOGICAL NUMBER CiF TAPE CONTAINING ELEMENT INPUT DATA FOR 
THIS SUBROUTINE. IF NUTEL - 5 t DATA IS READ FROM CARDS. 

NJ = NUMPER CF JOINTS OR ROWS IN MATRICES ( XYZ ) , (JDOF), (EUL). 

NUTMX = LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 
MASS MATRICES AND IVECS ARE OUTPUT. 

NUTMX MAY EE ZERO IF MASS MATRIX IS NOT FORMED. 

USFS FORTRAN PEAD, WRITE. 

NUTKX = LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 

STIFFNESS MATRICES (SAME AS GLOBAL LOADS TRANSFORMATION 
MATPICFS) AND IVECS ARE OUTPUT. 

NUTKX MAY PE ZERO IF STIFFNESS MATRIX IS NOT FORMED. 

USES FORTRAN READ, WRITE. 

NUTLT -• LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT LOCAL 



on on ooooooooooooooo on n n n r n 


AXIAL 


3/ 4 


LOAD TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. 

NUTLT MAY PE ZERO IF LOAD TRANSFORMATIONS ARE NOT FORMED. 
USFS FORTRAN READ, WRITE. 

NUTST = LOGICAL NUMEER OF UTILITY TAPE ON WHICH ELEMENT 

STRESS TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. 

NUTST MAY BE ZERO IF STRESS TRANSFORMATIONS ARE NOT FORMED. 
USES FORTRAN READ, WRITE. 

W = MATRIX WORK SPACE. MIN SIZE(6,6). 

T = MATRIX WORK SPACE. MIN SIZE(6,6). 

S = MATRIX WORK SPACE. MIN SIZE(6,6). 

KX = ROW DIMENSION OF XYZ IN CALLING PROGRAM. 

KJ = ROW DIMENSION OF JDOF IN CALLING PROGRAM. 

KF = ROW DIMENSION OF FUL IN CALLING PROGRAM. 

KW = ROW DIMENSION OF W, T, AND S IN CALLING PROGRAM. MXN=6. 

NEPROR EXPLANATION 

1 = JOINT NUMBER GREATER THAN NUMBER C c JOINTS. 

2 = MASS MATRIX FORMED, NUTMX . L E • ZERO. 

3 = STIFFNESS MATRIX FORMED, NUTKX .LE. ZERO. 

4 = LT MATRIX FORMED, NUTLT .LE. ZERO. 

5 = ST MATRIX FORMED, NUTST .LE. ZERO. 

1001 FORMAT (4<A6,4XJ) 

1002 FORMAT (2 (5X , E 10 .0 ) ) .. 

1003 FORMAT (3IE,4E10.0) 

2001 FORMAT (//46X 29H INPUT DATA FOR AXIAL ELEMENTS l 

2002 FORMAT (//40X 41HINPUT DATA FOR AXIAL ELEMENTS (CONTINUED)) 

2003 FORMAT (/ 16X7HMASS = A6, 9X7HSTIF = A6 , 9X13HL0AD TRANS = A6, 

* 6X1 5HSTRE SS TRANS = A6, 

* / 18X4HR0 = E10.3, 9X3HE = E10.3, 

* A/16X7HELEMENT 13X7HJ0INT 1 13X7HJ0INT 2 15X4HAREA 

* 16X4HAREA / 16X6H NUMEER 5 5X7HJ0INT 1 13X7HJGINT 2 /) 

2004 FORMAT (IX 3120, l^X E1G.3, 10X E10.3) 

READ AND WRITE FINITE ELEMENT DATA. 

NLINE = 0 
CALL PAGEHD 
WRITE (NOT, 2001) 

READ (NUT EL, 1 001 ) NAM EM ,NAMEK ,NAMELT,NAMEST 
READ (NUTEL»1 002) RO,E 

WRITE (NOT, 2003) NAME M, NAME K, NAME LT ,NAMEST,RO, E 
20 READ (NUTEL, 1 003 I NFL ,J1, J2,Al,A2 
IF (J1 .LE. 0) RETURN 
NLINE = NLINE + 1 
IF (NLINE .LE. 42) GO TO 30 
CALL PAGEHD 
WRITE (NOT, 2002) 

WRITE (NOT, 2003) NAMEM,NAME K,NAMELT LAMEST, RO ,E 
NLINE = 0 

30 WRITE (NOT, 2004) NEL, J1 , J2, A1 , A2 

NERR0R=1 

IF (J1 .GT. NJ .OR. J 2 .GT. NJ ) GO TO 999 

FORM FINITE ELEMENT COORDINATE LOCATIONS, EULER ANGLES REVADD IVEC. 
DO 42 1=1,3 
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CJCI.n = XYZfJlvI) 
CJ(I,2) = XYZ ( J2 * I ) 
EJ(I,1 ) = EUL(J1 ,1) 

E J ( I , 2 ) = FUL t J2 , I ) 
IV1 (I ) = JDC-F ( J 1 1 1 ) 

I VI (1+3 ) = JDOr- ( J2, I) 


c 

c 

c 


FORM MASS MATRIX (W). 

IF (NAMEM .EQ. 6H .CR. NAMEM .EO. 6HN0MASS) GO TG I1G 

CALL MAS1A (Cj,E J,A1 ,A2,RO, NAMEM, W,KCJ,KCJ,KW) 


IF (NUTMX .LE « 0) GO TO 999 

WRITE (NUTMX) NAMFM, N EL ,NPW ,NRW, NAMEL v ( IBLNK, 1=1 , 5) , 

* < CW(I,J ),I-1, NRW),J=1,NRW), (I VI (I), 1=1, NRW) 


NERRCR=2 


FORM STIFFNESS MATRIX (W), LOCAL LOAD TRANSFORMATION MATRIX (T), 
STRESS TRANSFORMATION MATRIX (S). 

110 IF (NAMEK .EQ. 6H .OP. NAMEK .EQ. 6HN0STIF) GO TO 20 

CALL STf-1 A ( C J » E J , A1 , A2 , F f NAMEK. »NAMEST ,W ,T ,S tNRST, 

* KCJ,KCJ,KW,KW»KW) 


I » W4\I\U * V" 


IF (NUTKX ,LE. 0) GO TO 999 
WRITE (NUTKX) NAMEK ,N EL ,NRK ,NRW ,NAMEL , ( IBLNK, I =1, 5 ) , 

* ( (W(I,J ), 1=1, NRW) ,J=1 ,NRW ) , ( IV1 (I ), 1=1, NRW) 

IF (NAMELT .EC. 6H .OR. NAMELT .EQ. 6HN0L0AD) GO TO 115 

IF (NUTLT .LE. C) GO TO 999 NERROR-4 

WRITE (NUTLT) NAMEL T, NE L, NRLT ,NRW,N AMEL , ( IBLNK , 1= 1, 5 ) , 

* ((T(I,J),I=1,NRLT),J = 1,NRW),(IV1(I) ,1-1, NRW) 

115 IF (NAMEST .EQ. 6H .OR. NAMEST .EQ. 6HN0STRS) GO TO 20 

IF (NUT ST .LE. 0) GO TO 999 NERR0R»5 

WRITE ( NUTST ) NAMEST, NE L , NRST , NRW , NAMEL ,( IBLNK , 1=1 , 5 ) , 

* (<S(I,J ),I=1,NRST),J=1,NRW),(IV1(I| ,I=l,Nf*W) 

GO TO 20 


C 


999 CALL ZZBOME (6HAXIAL ,NERROR) 
END 
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SUBROUTINE P1A1 (Rl -,Z»KZ) 

DIMENSION ZlKZ-1 S 

r 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C EUCKLING MATRIX 

C FOR AN AXIAL ROD ELEMENT WITH UNRESTRAINED F UNDARIE S. 

C BUCKLING MATRIX IS EASED ON UNIT AXIAL LOAD. 

C BUCKLING MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE POD TO lIF IN THE X-Z PLANE 
C WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

C LOCAL COORDINATE ORDER IS 
C DZ 1 »DZ2 

C WHFRE DZ IS TRAN S LA ‘ ION. 

C DEVELOPED BY RL WGHLEN. AUGUST 1973. 

C LAST REVISION SY RL WOHLEN. SEPTEMBER 19??. 

C 

C SUBROUTINE ARGUMENTS 

C RL = INPUT POD LENGTH. 

C 7 = OUTPUT EUCKLING MATRIX. SIZE (2*21. 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN-2. 

C 

C = 1 ./RL 
Z(ltl) ^ C 
Z ( 1 * 2 ) ~—C 

z(2,n =— c 

Z { 2 *2 ) = C 

RETURN 

END 
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SUBROUTINE B1A2 (RL,Z,KZ) 

DIMENSION Z(KZf)) 

C 

C SUER CUT ±NE TO CALCULATE FINITE ELEMENT.., 

C BUCKLING MATRIX 

C FOR A BEAM ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C BUCKLING MATRIX IS BASED ON UNIT AXIAL LOAD. 

C BUCKLING MATRIX IS IN LG CAL COORDINATE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE BEAM TO LIE IN THE X-l PLANE 

C WITH JOINT I AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

C LOCAL COORDINATE ORDER IS 
C DZ 1 ,DZ2 »TY1 ,TY2 

C WHERE DZ IS TRANSLATION AND TY IS ROTATION. 

C DEVELOPED BY RL WCHLEN. AUGUST 1973. 

C LAST REVISION BY RL WOHLEN. SEPTEMBER 1973. 

C 

C SUBROUTINE ARGUMENTS 

C RL = INPUT ROD LENGTH. 

C Z = OUTPUT BUCKLING MATRIX. SIZE(4,4). 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=4. 

C 

Cl = 6./(5.*RL) 

C2 = .1 

C3 = 12.*PL)/15. 

2(1,1) = Cl 
Z ( 1 ,2 ) =-C l 
Z(l,3) =-C2 
Z(l,4) =— C2 
Z C 2 ,2 ) = Cl 
Z ( 2 ,3 ) = C2 
2(2,4) = C2 
Z ( 2 , 3 ) = C3 
2(3,4) =-RL/30. 

Z ( 4,4) = C3 
DO 10 1=1,4 
DG 10 J=I,^ 

10 Z(J,I) = Z(I,J) 

C 

RETURN 

END 
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SUBROU. INF BAR ( XYZt JDOF -EUL , NUTEL ,N J» 

* NU TM X, NUTKX »NUTRX, NUTLT, NUTST, 

* W,T,S,KX,KJ,KE,KWI 

DIMENSION XYZ (KX, 1), JDOF i KJ, 1 ), EUL (KE,1 I tW (KW, 1 ) ,T(KW ,1 ) ,Sf KW ,1 1 
DIMENSION CJ( 3*3 ) f EJ 13,2 ), 1V1(12>, TR(12,12), TD(24,24) 

DIMENSION KODEK (4 ) , KODEB (2 ) « IFPIN(4), IV2I4I 
DATA NAMFL / 6HEAP / 

DATA N«?W,NRLT/12,12/, IBLNK/6H /, KCJ/3/, KTR/l2/,KTD/24/ 

DATA NIT ,NCT/5 ,6/ 

C 

C SUBROUTINF TO CALCULATE CON OPTION) FINITE ELEMENT ... 

C MASS MATRICFS AND IVECS (ON NUTMX). 

C STIFFNESS MATRICES (SAME AS GLOBAL LOAD TRANSFORMATION MATRICES) 

C AND IVECS (ON NUTKX )* 

C UNIT LOAD BUCKLING MATRICES AND IVECS (ON NUTBX), 

C LOCAL LOAD TRANSFORMATION MATRICES AND IVECS (ON NUTLT) * 

C STPFSS TRANSFORMATION MATRICES AND IVECS (ON NUTST) 

C FOR COMBINED AXI AL— TORSI ON-BENDING EAR ELEMENTS. 

C MASS, STIFFNFSS, BUCKLING MATRICES ARE IN GLOBAL COORDINATE 
C DIRECTIONS. 

C GLOBAL COORDINATE ORDER. FOR FACH ELEMENT IS 
C CU,V,W,P,Q,R> JOINT 1* THEN JOINT 2 

C WHERE U,V 7 W ARE TRANSLATIONS AND P,Q,R APE ROTATIONS. 

C IVEC GIVES ELEMENT DOF INTO GLOBAL DOF. EXAMPLES... 

C IVEC (6 ) =83*- PLACES ELEMENT DO c fc INTO GLOBAL DOF 834. 

C IVEC ( 3 ) =0 OMITS ELEMENT DCF 3 FROM GLOBAL DOF. THIS CONSTRAINS 

C ELEMENT DOF 3 TO ZERO MOTION. 

C GLOBAL LOAD TRANSFORMATION MATRIX RELATES LOAOS AT BAR ENDS IN GLOBAL 
C COORDINATE DIRECTIONS TO DEFLECTIONS IN THE GLOBAL COORDINATE 
C DIRECTIONS. 

C ROW ORDF° IN GLOBAL LOAD TRANSFORMATION MATRIX IS 
C (PU,PV,PW,MP,MO»MR) JOINT !, THEN JOINT 2. 

C WHERE P IS rf-RCE AND M IS MOMENT. 

C LOCAL LOAD T D ANSFOPMAT ION MATRIX RELATES LOADS AT BAP ENDS IN LOCAL 
C COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE DIRECTIONS. 

C ROW ORDER IN LOCAL LOAD TRANSFORMATION MATRIX IS 
C PX1,PX2,MX1,MX2,FY1,PY2,MZ1,MZ2,PZ1,PZ2,MY1,MY2 

C WHERE P IS FORCE AND MIS MOMENT. 

C STRESS TRANSFORMATION MATRIX RELATES STRESS AT BAR ENDS IN LOCAL 
C COORDINATE SYSTFM TO DEFLECTIONS IN THE GLOBAL COORDINATE DIRECTIONS. 

C ROW ORDER IN STRESS TRANSFORMATION MATRIX IS 
C PXI/AI ,PX2/A2» MX1*R1/TJ1,MX2*R2/TJ2, 

C PY1/A1 , PY2/A? »MZ 1*CY1 /B IZ 1 » MZ2* r V ?/BI Z2 , 

C PZ1/A1 ,PZ2/A2,MYl*CZI/ElYlyMY2*CZ2/BIY2 

C WHERE P IS FORCE AND M IS MOMENT. 

C DATA ARRANGEMENT ON NUTMX, NUTKX, NUTtX, NUTLT, NUTST FOR EACH 
C FINITE ELEMENT IS ( W=M ,X ,B ,LT , ST ) 

C WRITE ( NUTWX ) NAM E W ,N FL ,NR , NC ,NAMEL , ( IB LNK, 1=1 , 5) , 

C ( (W(I,J),1=1,NR),J=1,NC),(IVEC(I) , I =1 ,NC ) 

C CALLS FORMA SUBROUTINES EUC1E, MAS IE. , PAGE HD, STF18, ZZBOME. 

C DEVELOPED BY WA EENFIELD, CS BODLEY, RL WOHLEN. FEBRUARY 1973. 
r LAST REVISION BY RL WOHLEN. APRIL 1976. 

C 

C ♦********¥**£****#********* ***** ****************************** ************ 

C INPUT DATA READ IN THIS SUBROUTINE FROM NUTEL. IF NUTEL = 5, DATA IS 
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READ FROM CARDS. 

NAMEM,NAMFK t NAMLLT,NAKEST,NAMEe, 

(KODEK Cl ),I=1*4),(K0DEB(I)»1=1»2) 

RO,E,C 

20 NEL,J1, J2,JREF,A1,PI1 ,TJ1 ,BIZ 1 ,BIY1 ,SF* 
IFPY1,IFPZ1,IFPY2,IFP Z2,XFTAPR 
IF (J! .FC. 0) RETURN 

IF (IFTAPR .EO. 1HT) A2,PI2 ,TJ2,B IZ 2,81 Y2 
30 IF (LAMEST .EC. 6h .OR. NAMEST .EO. 

R1 ,CY1,CZ’,R2,CY2,CZ2 
GO TO 20 


FORMAT (5(A6*4X) , 
A1,A1»A?,A2,4X,A2,A2) 
FORMAT (3(5X,E1C)) 

FORMAT (4I5,5E10,E5 ,5A1 ) 

FORMAT (2CX,5E10) 
6HN0STRS1 GO TO 20 
FORMAT (20X,6E 10) 


INPUT DATA REQUIREMENTS 

AXIAL 

TORSION 

BENDING 

ALCNG 

ABOUT 

AEOUT 

LOCAL X 

LOCAL X 

LOCAL Z 


BENDING 
ABOUT 
LOCAL Y 


MASS 

STIF, LOAD TRANS 
EUCKLING 
STRESS TRANS 



PI»PO A,PO A,RO 

TJ,G EIZ,A*SF,E,G BIY,A,SF,E,G 

NONE NONE NONE 

STIF+R STIF^CY STIF+CZ 

FOR NO SHEAR DEFC’RMATIGW IN BENDING, SET ANY OF A (NOT IF AXIAL USED), 
SF, OP G (NOT IF TORSION IS USED) TO ZERO. IF BENDING STRESS 
TRANSFORMATION IS WANTED, A MUST NOT BE ZERO. 


DEFINITION OF INPUT VARIABLES. 

NAMEM - TYPE OF MASS MATRIX WANTED. 

= Ml, DIAGONAL LUMPED. 

= M2, CONSISTENT. 

= th Cfc 6FN0MASS, NO MASS MATRIX CALCULATED. 

NAMEK = TYPE OF STIFFNESS MATE 1 X WANTED. 

= Kl* CONSTANT FORCE FOR AXIAL, CONSTANT TORQUE FOR TORSION, 
CONSTANT SHEAR AND LINEAR MOMENT FOR BENDING. 

= 6H OF fcHNOSTIF, NO STIFFNESS MATRIX CALCULATED. 

NAMELT = IDENTIFICATION NAME FIR LOAD TRANSFORMATION MATRICES. 

= 6H CR 6FN0LCAD, NO LOAD TRANSFORMATIONS CALCULATED. 

NAMEST - IDCNTIFICATION NAME FOR STRESS TRANSFORMATION MATRICES. 

= 6H OR fcHNOSTRS, NO STRESS TRANSFORMATIONS CALCULATED. 

NAMED = TYPF OF EUCKLING MATRIX WANTED. 

- FI, AXIAL RCO. 

= F 2 , BEAN . 

= fcp OR 6 KNTPUCK, NO BUCKLING MATRIX CALCULATED. 

KODEK = OPTION COCF FOP AXIAL, TORSION, BENDING Z, BENDING Y LOCAL 
STIFFNESS. IF BLANK, 'LL FOUR ARF CALCULATED. SIZE (4). 

KODEK (1 J-A , LOCAL STIF MATPIX IS CALCULATED FOR AXIAL 
(ALONG LOCAL X-AXIS). 

KODEK (2 )=T , LOCAL STIF MATRIX IS CALCULATED FOR TORSION 
(AFOOT LOCAL X-AXIS? «_ 

KDDEK (3 )=BZ» LGCAl STIF MATRIX IS CALCULATED FOR BENDING 
(AFCUT LOCAL Z-AXISI. ^ 

KODFK(A-)=BY, LOCAL STIF MATRIX IS CALCULATED FOR BENDING 
(ABOUT LOCAL Y-AXIS). 

KODEB = OPTION CODT FOR BUCKLING IN LOCAL Y CR Z DIRECTION. 

IF FLANK, PCTH ARF CALCULATED. SIZF(2). 

KODEB (I )=BY, LOCAL EUCKLING MATRIX IS CALCULATED FOR 
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DEFLECTION IN LOCAL Y DIRECTION. 

KODEE (2 )=E Z , LOCAL BUCKLING MATRIX IS CALCULATED FOR 
DEFLECTION IN LOCAL Z DIRECTION. 

RO = MASS DENSITY. 

E = YOUNGS MODULUS OF ELASTICITY. 

G = SHEAR MODULUS OF ELASTICITY. 

NEL = FINITE ELEMENT NUMBER. FOR REFERENCE ONLY, NOT USED IN 
CALCULATIONS. WRITTEN ON NUTMX, ETC. 

Jl = JOINT NUMBER AT BAR END 1. LOCAL X-AXIS ORIGINATES AT Jl. 

J2 = JOINT NUMBER AT BAR END 2. LOCAL X-AXIS GOES FROM Jl TO J2. 

JREF = REFERENCE POINT. LOCAL Z-AXIS IS DEFINED BY VECTOR (J1,J2> 

CROSSED INTO VECTOP tJl,JREF). LOCAL Y-AXIS LIES IN XY PLANE 
DEFINED EY J1,J2,JPEF. 

A1 ~ CROSS-SECTION AREA AT BAR END 1. 

A2 = SAME AS Al AT BAR End p. 

PI1 = CROSS-SECTION POLAR AREA MOMENT OF INERTIA FOR MASS 
CALCULATIONS AT EAR FND 1. 

PI2 = SAME AS Pll AT EAR END 2. 

TJ1 = CROSS-SECTION SAINT VENANTS TORSION CONSTANT (Jl IN JG FOR 
TORSION STIFFNESS AT EAR END 1. 

TJ2 = SAME AS TJl AT BAP END 2. 

BIZI = CROSS— SECTION AREA MOMENT OF INERTIA ABOUT LOCAL Z-AXIS 
(FC* BENDING) AT FAR END 1. 

BIZ2 = SAMF AS FIZ1 AT EAR END 2. 

BIYI = CROSS— SECTION AREA MOMENT OF INERTIA ABOUT LOCAL Y-AXIS 
(FCR FENCING! AT PAP END 1. 

6IY2 ~ SAME AS BIYI AT EAR FND 2. 

SF = SHAPE FACTOR (K I F"P SHEAR IN KAG. 

USE SF=0.0 FOR NO SHEAR DEFORMATION IN BENDING. 

SF= 1.0 FOP A SOLID CIRCULAR CYLJNDFR. 

SF= .5 FOP A THIN WALLED CIRCULAR CYLINDcR. 

IFPY 3 = PIN JOINT OPTION FOR LOCAL COORDINATE THETA Y AT BAR END 1. 

= 1H , MOMENT JOINT. 

= 1 HP, PIN JOINT. 

IFPY2 = SAME AS IFPY1 AT EAR END 2. 

IFPZ1 = PIN JOINT OPTION FOR LOCAL COORDINATE THETA Z AT BAR END 1. 

= 2H , MOMENT JOINT. 

- InP, PIN JCINT. 

IFPZ2 = SAME AS IFP21 AT FAR FND 2. 

IFTAPR = OPTION FOP TAPERED EAR. 

= 1H , CONSTANT SECTION PROPERTIES. 

= 1 HT, LINEAR TAPER SECTION PROPERTIES. 

RI = DISTANCE FROM LOCAL X-AXIS TO CUTER FIBER FOR TORSION 

STRESS CALCULATION AT EAR END 1. 

R2 = SAMF AS RI AT B A c END 2. 

CY1 .= 0 1 STANCE FROM LOCAL XZ PLANF TO OUTER FIBER FOR BENDING 

STRESS CALCULATION AT BAR ENU 1. LOCAL Y DIRECTION. 

CY2 = SAMF AS CY1 AT BAR END 2. 

CZ1 = DISTANCE FROM LOCAL *Y PLANE To OUTER FIBER FOR BENDING 

STRESS CALCULATION a I BAP END 1. LOCAL Z DIRECTION. 

CZ2 = SAME AS C7.1 AT BAR END 2. 

EXPLANATION OF INPUT FORMATS. NUMBER INDICATES CARD COLUMNS USED. 

I = INTEGER DATA* RIGHT ADJUSTED. 

E = DECIMAL POINT DATA, ANYWHERE IN FIELD. EXPONENT RIGHT ADJUSTED. 
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X = CARD COLUMNS SKIPPED. 

SUBROUTINE ARGUMENTS CALL I.NPUT) 

XYZ a MATRIX OF JOT NT GLOBAL X,Y,Z LOCATIONS. ROWS CORRESPOND 
TC JOINT NUMBERS. COLUMNS 1,2*3 CORRESPOND TO THE JOINT 
X f Y *Z LOCATIONS RES ECTIVELY. STZE(NJ,3). 

JDOF = MATRIX OF JOINT Gr'U' AL DEGREES OF FREEDOM. ROWS CORRESPOND 
TO JOINT NUMBER J.U PLUMNS 1*2,3 CORRESPOND TO ThE JOINT 
TRANSLATION DOFS A/D COLUMNS 4,5,6 CORRESPOND TO THE JOINT 
ROTATION OOFS. SXZ;(NJ,6I. 

EUL = MATRIX OF JOIAT rttLER ANGLES (DEGREES!. ROWS CORRESPOND 
' TO JOINT NUMBERS.: COLUMNS 1,2,3 CORRESPOND TO THE 
GLOBAL X ,Y,Z PERMUTATION. SIZE(NJ,3). 

NUTEL = LOGICAL NUMBER OF' TAPE CONTAINING ELEMENT INPUT DATA FOR 
THIS SUBROUTINE. IF NUTEL = 5, DATA IS READ FROM CARDS. 

NJ = NUMBER OF JOINTS tjf ROWS IN MATRICES (XYZf, (JDOF), (EUL). 

NUTMX =■ LOGICAL NUMBER OF l TILITY TAPE ON WHICH ELEMENT 
MASS HA TP ICE'S AND IVECS ARE OUTPUT. 

NUTMX MAY BE ZERO IF MASS MATRIX IS NOT FORMED. 

USES FORTRAN READ, WRITE. 

NUTKX = LOGICAL NUMBER 1 OF UTILITY TAPE ON WHICH ELEMENT 

STIFFNESS MATRICES ISAMF AS GLOBAL LOADS TRANSFORMATION 
MATRICES) AND I VECS ARE OUTPUT. 

NUTKX MAY EE ZERO IF STIFFNESS MATRIX IS NOT FORMED. 

USES FORTRAN REA V >, WRITE. 

NUTBX = LOGICAL NUMBER Ol : UTILITY TAPE ON WHICH ELEMENT UNIT LOAD 
BUCKLING MATRICES AND I VECS ARE OUTPUT. 

NUTBX MAY BE ZERO IF BUCKLING MATRIX IS NOT FORMED. 

USES FORTRAN REAL* WRITE. 

NUTLT =• LOGICAL NUMBER Of UTILITY TAPE ON WHICH ELEMENT LOCAL 
LOAD TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. 

NUTLT MAY EE ZERO IF LOAD T' *NS FORM AT IONS ARE NOT FORMED. 
USES FORTRAN READ, WRITE. 

NUTST = LOGICAL NUMBER CF UTILITY TAPF ON WHICH ELEMENT 

STRESS TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. 

NUTST MAY EF ZERO IF STRESS TRANSFORMATIONS ARE NOT FORMED. 
USES FORTRAN READ, WRITE. 

W = MATRIX WORK SPACE. MIN SIZFf12*12). 

T = MATRIX WORK SPACE. MIN $IZE(12?12:. 

S = MATRIX WORK SPACE. MIN SIZE (1 0 >T2 ) . 

KX = ROW DIMENSION CF XV 2 I,N CABLING PROGRAM. 

KJ = ROW DIMENSION OF JDOF IN C, uLIN'G PROGRAM. 

KE = ROW DIMENSION OF EUL IN -aLLING PROGRAM. 

KW = ROW DIMENSION OF W, T, AND S IN CALLING PROGRAM. HIN=12. 

NERRCF EXPLANATION 

1 = JOINT NUMBER GREATER T AN NUMBER OF JOINTS. 

2 = STIFFNESS MATRIX FOR ' ... p, NUTKX .LE. ZERO. 

3 = LOAD TRANSFORMATION MATRIX FORMED, NUTLT .LE. ZERO. 

4 = STRESS TRANSFORMATION MATRIX FORMED, ''iUTST .LE. ZERO. 

5 = MASS MATRIX FOkMFO, NUTMX -L c . ZFRO. 

6 = BUCKLING MATRIX FORMED, NH-’BX .LE. ZERO. 

1001 FORMAT <5(A6,4X),Al*Al m2 ,A2 , 4X A2,A2) 
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1002 FORMAT (3(5X,E10.G) ) 

1003 FORMAT (* 15 ,5F10.0, E5 .0 ,5 Al > 

1004 FORMAT (20X,6E1C.0) 

2001 FORMAT I//46X 27HINPUT DATA FOR BAR ELEMENTS) 

2002 FORMAT C//40X 39HINPUT DATA FOR BAR ELEMENTS (CONTINUED)) 

2003 FORMAT ( 45X,8HK0DEK = A1*A1*A2,A2* 4X 8HK0DEB * A2,A2* 

c * / 10X7HMASS = A6* 6X7HSTIF = A6» 6X13HL0AD TRANS = Afe, 

* 3X15HSTRESS TRANS = A6- 3X1 1HBUCKLING = A6* 

* / 12X4HRC = E10.3, 6X3HE = E10.3, 80X7HI III* 

* / 32X3FG = E10.3, 80X7HF F F F, 

* / 125X7HP P P P, 

* / 1X7HELEMENT 2X5HJ0I NT 2X5HJ0INT 3X3WEF 5X4HAREA 

'* 7X5HP0LAP. 5X7HTORSION 3X9H2 BENDING 2X9HY BENDING 

* 2X5HSHEAP 3X6HSTF ESS 5X6H STRESS 5X6HSTRESS 3X7HY Z Y Z 

* / 1X6HNUMBER 5X1H1 6X1H2 4X5H POINT 

* I4X 7HINERTIA 5X5HC0NST 5X7HINEPTIA 

* 4X7HINERTIA 3XBHF ACTOR 5X1HR 9X2HCY 9X2HCZ 5X7H1 1 2 2/1 

2004 FORMAT (IX 15, 18, 217, IX 5E11.3* F7.3, 3E11.3, 4(1XA1)) 

2005 FORMAT (29X,5E11 .3,7X ,3E1 1.3) 

C 

C READ AND WRITE FINITE ELEMENT DATA._ 

R 1 = 0.0 
CY1 = 0.0 
CZ1 - 0.0 
NLINE = 0 
CALL PAGEFD 
WRITE (N0T,2001) 

READ (NUTEL,1001) NAMEM,NAMEK,NAMELT,NAMEST,NANEB, 

* (KODEK (I) ,1=1,4), (KOCEBd )*I=1*2) 

RFAD (NUTEL, 1002 ) RO, E,G 

WRITE (N0T,2GC3) (KCDEK(I ) , 1=1 ,4) , ( KODEB ( I ) ,1=1,2), 

* NAME M* NAME K* NAME LT *NAMEST*NAMEB»RO*F *G 
20 READ (NUTEL,1003) NEL ,J1, J2*JREF,Al,Pll,TJl,BI2l,BIYI,SF, 

* I FP IN, I FT APR 
IF (J1 .LE. 0) RETURN 

IF (IFTAPR .EQ. 1HT) READ (NUTEL,1004) A2,P12 ,TJ2 ,BIZ2,BIY2 
IF (NAMFST .TO. 6h .OR. NAMEST .EO. 6HN0STRS) GO TO 25 

READ (NUTFL, 100*7 R1,CY1,C21,R2,CY2,C22 
25 NLINE = NLINE ♦ ) 

IF (IFTAPR .EC. 1HTI NL INF=NL INE+1 
IF (NLINE .LE. 42) GO TO 30 
CALL PACE PC' 

WRITE (NOT, 2002) 

WRITE *NCT*?G03) (KCD EK ( 1 ) , 1=1 ,4) * (KOOEE ( I) *1=1*2 ). 

* NAMEM,NAMEK,NAMELT,NAMEST,NAKEB,RO,E,G 
NLINE = C 

30 WRITE (NO! ,2004) NE L, JJ , J2 , JR EF ,A1 ,PI 1,TJ1, BIZI *B IY1 ,SF , 

* R1 ,CYj »C21 ,IFP1N 

NERR0R=1 

IF (J1 .&T. NJ .OR. J 2 .GT. NJ .OR. JREF .GT. NJ) GO TO 999 
IF (IFTAPR .FO. 1HT) WRITE (N0T,2005» A?, PI 2, TJ2* B1Z2 ,& IY2 »R2 * 

* CY2,CZ2 

FORM FINITE FLEMENT COORDINATE LOCATIONS, EULER ANGLES* REVADD IVEC. 
DO 42 1=1,3 
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CJCI.l) = XY2CJ1.II 
CJCI.2) = XY2t,»2,n 
CJCI.3) = XY2CJREF V I} 

EJCI.l) = EULCJltl) 

42 FJCI.2) * EUL C J2 . X 1 
DO 44 1*1 .6 
IV1CIJ = JDOFCJ1.I) 

44 IV1CI+61 * JDOFCJ2.I) 

FORM DATA FOP UNIFORM ELEMENT. 

IF CIFTAPR .EC. 1HT) GO TO 50 

A2 = A1 

PI2 = PI1 

TJ2 = TJ1 

BI22 = PI21 

BIY2 = 6IY1 

R2 = R1 

CY? = CY1 

C22 = CZ1 

FORM PINING IVEC. 

50 NPIN = 0 


DO 

55 

1=1, 

L. 





IF 

CIFPINC 

I) , 

»NE. 

1 HP ) 

GO 

TO 

NPIN = 

NPIN ♦ 

I 




IF 

Cl 

.EQ. 

1) 

IV2 

(NPIN) 


11 

IF 

(I 

.EO. 

21 

IV2 

CNPIN) 

— 

7 

IF 

Cl 

.EO. 

3) 

IV2 

CNPIN) 


12 

IF 

Cl 

.EQ. 

4) 

1V2 

CNPIN) 


8 

55 CONTINUE 







FORM STIFFNESS MATRIX CM), LOCAL LOAD TRANSFORMATION MATRIX CT1* 
STRESS TRANSFORMATION MATRIX CS). 

100 IF (NAMFK .EO. 6H .OP. NAME* .EQ. 6HN0STIF! GO TO 110 

CALL STF1B CCJ.EJ.KODEK.Al .A2.T II .TJ2.BIZ1 »B 122, BI Y1 ,6 IY2 *R1 ,R2, 

* CY1 ,CY2, CZ 1 «CZ2 , ST »E .G.NAMEK.NAMEST.K.T.S.NRST , 

* KCJ.KCJ.KW.KW. KWJ 

PIN STIFFNESS MATRIX. 

IF (NPIN .FQ. C) GO TO 105 
CALL DCPS1B (CJ.EJ.H, KCJ.KCJ.KW) 

CALL TRANS CW.TP ,NRW ,NRW .KW.KTR 1 
CALL MULTA CT.TR.NRLT, NRW, NRW, KW.KTR) 

CALL SPED3 C T , IV2 »T, TD .NRW.NPIN ,1 , KW > 

CALL MULTA ( TD, W ,NPW ,NRW ,NRW »KTD .KW) 

CALL MULTA ( TP ,TD,NR W, NRW, NR W.KTR , KTD) 

CALL MULTA ( T.W ,NRLT ,NRW ,NRW .KW.KW 1 
CALL ATXRA1 CW.T.NRW.NRW.KW.KWl 

IF (NAME ST .EC. 6H .OR. NAMEST .EO. 6HN0STRS) GC TO 105 

CALL MULTA C S.TR »NRST .NRW.NR W.KW ,K TR ) 

105 N£RR0R=2 

IF (NUTKX .LE. 0) GO TO 999 

WRITE (NUTKX) NaMFK ,NEL, NRW »NRW,NAMEL ,( IBLNK, 1=1, 5» , 

* ((WCI,J),I=1,NRW)»J=1,NRW),CIV1(I), I-l.NRW) 

IF CNAKELT .EQ. fcH .OR. NAMELT .EC. 6HN0L0AD) GO TO 115 
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NERR0R=3 

IF CNUTLT .LE. 0) GO TO 999 

WRITE (NUTLT) NAMELT,NFL,NRLT,NRW,NAMEL , (IBLNK, 1=1,5) , 

* ((T(ltJ)*l=lf!*LT) t J=l,NRW),(IVim ,I=1,NRW) 

115 IF CNAMEST .EO. 6H .OR. NAHEST »EC. 6HN0STRS) GO TO 110 

NERR0R=4 

IF (NUTST .LE. C) GO TG 999 

WRITE {NUTST) NAMEST»NEL,NRST ,NRW,NAMEL ,(IBLNK,I=1,5), 

* ((S<I,J),1=1,NRST),J=1,NRW),«IV1CI) ,1=1, NRW) 

FORM MASS MATRIX (W). 

110 IF (NAMEM .EO. 6H .O r *. NAMEM .EO. 6HN0MASS) GO TO 140 

CALL MAS1B 1C J,EJ, Alt A2,PI 1»PI2,R0,NAMEM,W,T,KCJ ,KCJ,KW,KW) 

PIN MASS MATRIX. 

- IF CNPIN .GT. 0) CALL BTABA (W,TR,NRW,NRW,KW,KTR ) 

NERR0R=5 

IF CNUTMX .LE. 0) GO TO 999 

WRITE CNUTMX) NAMEM ,NEL,NRW * NRW, NAMEL .( IBLNK, 1=1, 5) , 

* ( (W (I ,J ) , 1=1 » NRW) * J=1 1 NRW ) * (I VI Cl ), 1=1, NRW) 

FORM UNIT LOAD BUCKLING MATRIX (W). 

140 IF (NAMEB .EO. 6H .OR. NAMEB .EO. 6HN0BUCK) GO TO 20 

CALL BUC1B (CJ,E J,KO DEB, NAMEB, W,S ,KCJ,KCJ,KW,KW) 

NERROR-6 

IF CNUTPX .LE. 0) GO TO 999 , 

WRITE (NUTBX) NAMEB ,N EL , NRW , NRW, NAMEL ,( IBLNK, 1=1, 5) y 

* ( (W (I , J ), 1=1, NRW) ,J=1 ,NRW ),(IV1(I), 1=1, NRW) 

GO TO 2C- 

999 CaLL Z2B0MB (6HBAR ,NERROR) 

END 
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SUBROUTINE BUC1B fCJ,EJ, KODEB, NAMEB,Z,W,KCJ,KEJ,KZ,KWJ 
DIMENSION CJ(KCJ,1>, EJ(KEJ,1), KODEBC1 1* ZCKZ,1| t W(KW,U 

SUBROUTINE TO CALCULATE FINITE ELEMENT... 

BUCKLING MATRIX 

FOR A COMBINED AXIAL —TOR SION -BENDING BAR ELEMENT WITH UNRESTRAINED 
BOUNDARIES. 

BUCKLING MATRIX IS BASED ON UNIT AXIAL LOAD. 

BUCKLING MATRIX IS IN GLOBAL COORDINATE DIRECTIONS. 

GLOPAL COORDINATE ORDER IS 

(U,V,W,P,0,R) JOINT 1, THEN JOINT 2 
WHERE U,V,W ARE TRANSLATIONS AND P,0,R ARE ROTATIONS. 

EULER ANGLE CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

CALLS FORMA. SUBROUTINES B1A1 , BIA2* BTABA* DCOS1B, ZZBOMB. 

DEVELOPED BY RL WOHLEN. AUGUST 1973. 

LAST REVISION BY WA BENF1ELD. MARCH 1976. 

SUBROUTINE ARGUMENTS 

CJ = INPUT MATRIX OF GLOBAL X,Y,Z COORDINATES AT BAR JOINTS. 

ROWS 1,2,3 CORRESPOND TO X,Y,Z COORDINATES. 

COLS 1,2 CORRESPOND TO JOINTS 1,2. COL 3 CORRESPONDS 
TO REFERENCE POINT TO DEFINE LOCAL XY PLANE. SIZE(3,31. 
EJ = INPUT MATRIX OF EULER ANGLES CDEGREFS) AT BAR JOINTS. 

ROWS 1,2,3 CORRESPOND TO GLOPAL X,Y,Z PERMUTATION. 

COLS 1,2 CORRESPOND TO JOINTS 1,2. SIZE<3,2). 

KODEB = INPUT OPTION CODF FOR LOCAL Y, LOCAL Z BUCKLING. 

IF BLANK, BOTH ARE CALCULATED. SIZE(?I. 

KCDEB ( 1 ) =BY , LOCAL BUCKLING MATRIX IS CALCULATED 
FOR LOCAL Y DIRECTION. 

KCDE6( 2 ) =LZ, LOCAL eUCKLING MATRIX IS CALCULATED 
FOR LOCAL Z DIRECTION. 

NAME6 = INPUT TYPE OF BUCKLING MATRIX WANTED. 

=P1, AXIAL ROD. 

= 62 , e FAM. 

Z = OUTPUT BUCKLING MVTRIX. SIZFU2,12i. 

W = INPUT WORK SPACE MATRIX. SIZF(12,12). 

KCJ = INPUT ROW DIMENSION OF CJ IN CALLING PROGRAM. 

KEJ = INPUT ROW DIMENSION OF EJ IN CALLING PROGRAM. 

KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=12. 

KW = INPUT ROW DIMENSION OF W IN CALLING PROGRAM. MIN=12. 

NERRPF. EXPLANATION 

1 = DIMENSION SIZE EXCEEDED. 

2 = IMPROPERLY DEFINED NAMEB. 

NERROR=l 

IE (KZ .IT. 12 .OR. KW .LT. 12) GO TO 999 
DO 6 J*I,1? 

DO 5 1=1,12 
5 Z(I,J) = 0.0 

P-L = SQPT<«CJ(1,2)-CJ U ,1 >1**2 + (C J ( 2, ? )-C J ( 2 ,1) )**2 
* ♦ (CJ(3,2)-CJ(3,1) )**2 ) 

KODEBY = 1 
KODEPZ = 1 

IF (KODEB (1 ) .E0.2H .AND. KODEB ( 2 ) . EO .2H ) GO TO 10 
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IF (KODEetl) .NE. 2HBY) KODEBY = 0 
IF (KODEB(2) .NF. 2HBZ) KODEBZ = 0 
10 IF (NAMEE .EQ. 6HB1 ) GO TO ’10 

IF (NAMEB .EQ. 6HB2 I GO TO 120 

NERR0R=2 

GO TO 999 

110 IF (KODEBY .EO. 1) CALL BlAl (RL,Z ( 5,5 1 ,KZ) 

IF (KODEB Z .EC. 1) CALL BlAl (RL,Z ( 9,9) ,KZ) 

GO TO 300 

120 IF (KODEBY .EC. 1) CALL B1A2 (RL,Z ( 5,5) ,KZ) 

DO 125 J=7,8 
DO 125 1=5,6 
Z(1,J) =-Z(l,J> 

125 Z(J,I) =-Z(J,I> 

IF (KODEBZ .EQ. 1) CALL B1A2 (RL, 2(9,9) ,KZ) 

300 CALL DC0S1R (CJ, EJ,W, KC J,KEJ,KW) 

CALL BTABA (Z,W, 12,12, K2,KW) 

RETURN 

999 CALL ZZBOMB (6HEUC1B ,NERROR) 

END 
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SUBROUTINE DCOSIA (CJ ,EJ, 2 ,KC J,KEJ ,KZ ) 

DIMENSION CJ(KCJ,1), EJ(KEJ,1), Z(KZ,1) 

DIMENSION P(3 ) * T(3,3) 

SUBROUTINE TO CALCULATE FINITE ELEMENT... 

DIRECTION COSINE MATRIX 
FOR AN AXIAL ROD ELEMENT. 

THE DIRECTION COSINE MATRIX RELATES LOCAL COORDINATE DISPLACEMENTS 
TO GLOBAL COORDINATE DISPLACEMENTS. 

THE LOCAL COORDINATE SYSTEM ASSUMES THE POD TO LIE ALONG THE X AXIS 
WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

ROW ORDER (LOCAL COORDINATE ORDER) OF DIRECTION COSINE MATRIX IS 
DX1 »DX? 

WHERE DX IS TRANSLATION. 

COLUMN ORDER (GLOBAL COORDINATE ORDER) OF DIRECTION COSINE MATRIX IS 
(U,V,W) JOINT 1, THEN JOINT 2. 

WHERE U,V,W ARE TRANSLATIONS. 

EULER ANGLE CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

CALLS FORMA SUBROUTINES EULER ,MULTB,ZZBOMB . 

DEVELOPED BY RL KCHLEN. SEPTEMBER 1972. 

LAST REVISION BY WA BENFIELD. MARCH 1976. 

SUBROUTINE ARGUMENTS 

CJ = INPUT MATRIX OF GLOBAL X,Y,Z COORDINATES AT ROD JOINTS. 

ROWS 1,2,3 CORRESPOND TO X,Y,Z COORDINATES. 

COLS 1,2 CORRESPOND TO JOINTS 1,2. SIZE(3,2). 

EJ = INPUT MATRIX OF EULER ANGLES (DEGREES) AT ROD JOINTS. 

ROWS 1,2,3 CORRESPOND TO GLOBAL X,Y,Z PERMUTATION. 
CCLS 1,2 CORRESPOND TO JOINTS 1,2. SIZEC3,2). 

Z = OUTPUT DIRECTION COSINE MATRIX. SIZE(2,6). 

KCJ = INPUT ROW DIMENSION OF CJ IN CALLING PROGRAM. 

KEJ = INPUT ROW DIMENSION OF EJ IN CALLING PROGRAM. 

KZ - INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=2. 

NcP RPR EXPLANATION 
1 = DIMENSION SIZE LESS THAN 2. 

NERROR = 1 

IE (KZ .LT. 2) GO TO 999 
PX = CJ(1 ,?)-CJ(I ,1 ) 

PY = CJ ( 2 ,2 )— C J ( ? ,1 ) 

PZ = CJ (3 ,2 )— C J (3 ,1 ) 

PL = SORT (PX**2 + PY**2 + PZ**2) 

P ( 1 ) = PX/PL 
P (2 ) = PY/PL 
P ( 3 ) = PZ/PL 
DO 10 1=1,2 
DO 10 J = 1 ,6 
10 Z(I,J) = 0.0 

CALL EULFP ( E J ( 1 , 1 ) , T, 3) 

CALL MULTB (P,T, 1,3,3, 1,3) 

DO 22 J=1 ,3 
22 Z(1,J) = T ( 1 , J ) 

CALL EULFP ( E J ( 1 ,2 ) , T, 3 ) 

CALL MULTB (P,T, 1,3,3, 1,3) 
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DO 24 J=l,3 
24 Z (2* J+3) = 

RETURN 

999 CALL ZZtsOMB ( 6HDC0S1A »NERROR) 
END 
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C 

C 

C 

C 

r 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE- DCOSlri ( CJ ,E J , Z ♦ KC J *KE J , KZ ) 

DIMENSION CJ(KCJ»1), EJ IKEJ , 1 ) , Z (KZ, 1) 

DIMENSION W(3,3), T(3,3) 

SUBROUTINE TO CALCULATE FINITE ELEMENT... 

DIRECTION COSINE MATRIX 

FOP A COMBINED AXIAL-TORSION-BENDING BAR ELEMENT. 

THE DIRECTION COSINE MATRIX RELATFS LOCAL COORDINATE DISPLACEMENTS 
TO GLOBAL COORDINATE DISPLACEMENTS. 

THE LOCAL COORDINATE SYSTEM ASSUMES THE BAR TO LIE IN THE X-Y PLANE 
WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS, 
REFERENCE POINT 3 (TO DEFINE THE LOCAL X-Y PLANE) IS IN THE 
POSITIVE Y DIRECTION. 

ROW ORDER (LOCAL COORDINATE ORDER) OF DIRECTION COSINE MATRIX IS 
DX1 *DX2 , TX1 , TX2 , DY1 ,DY2 ,TZ1 ,TZ2, D2 1, DZ2, TY1 ,TY2 
WHERE DX,DY,DZ ARE TRANSLATIONS AND TX,TY,TZ ARE ROTATIONS. 

COLUMN ORDER (GLOBAL COORDINATE ORDER) OF DIRECTION COSINE MATRIX IS 
(U,V,W,P,Q,P) JOINT 1, THEN JOINT 2 
WHERE U,V,W ARE TRANSLATIONS AND P,Q,R ARE ROTATIONS. 

EULER ANGLE CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

CALLS FORMA SUBROUTINES EULER, MULTB ,ZZBOMB . 

DEVELOPED BY RL WOHLEN. FEBRUARY 1973. 

LAST REVISION BY WA BENFIELD. MARCH 1976. 


SUBROUTINE ARGUMENT S 

CJ = INPUT MATRIX OF GLOBAL X,Y,Z COORDINATES AT BAR JOINTS. 

ROWS 1,2,3 CORRESPOND TO X,Y,Z COORDINATES. 

COLS 1,2 CORRESPOND TO JOINTS 1,2. COL 3 CORRESPONDS 
TO REFERENCE POINT TO DEFINE LOCAL XY PLANE. SIZE(3,3) 
EJ = INPUT MATRIX OF EULER ANGLES (DFGPEES) AT BAR JOINTS. 

ROWS 1,2,3 CORRESPOND TO GLOBAL X,Y,Z PERMUTATION. 

COLS 1,2 CORRESPOND TO JOINTS 1,2. SIZE(3,2). 

Z = OUTPUT DIRECTION COSINE MATRIX. SIZE(12,12). 

KC J = I h'P UT ROW DIMENSION OF CJ IN CALLING PROGRAM. 

KEJ = INPUT ROW DIMENSION OF EJ IN CALLING PROGRAM. 

KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=12. 


NERROR EXPLANATION 
1 = DIMENSION SIZE LESS THAN 12. 


NERROR - I 

IF (KZ .LT. i2 ) GO TO 999 
PX = CJ(1 ,2)— CJ(1 ,1 ) 

PY = C J ( 2 ,2 J-C J (2,1 ) 

PZ = C J (3,2 )-C J( 3,1 ) 

PL = SORT (PX**’ + PY**2 +• PZ**2) 

RX = PY*(CJ(3»3)-CJ (3,1 )) - PZ* < C J ( 2, 3 ) -C J { 2, 1 )) 

PY = PZ*(CJ(1 ,3)-CJ(l ,1 )} - PX*(CJ(3,3)— CJ(3,1 )) 

RZ = PX*(CJ(2,?)-CJ (2 ,1 )) - PY*(CJ (1*3)— CJ(lfl)) 

RL = SORT (RX**2 + RY**2 + R2**2) 

OX = RY*PZ - RZ*PV 
OY = R Z *P X - PX*P2 
CZ = RX*PY - RY*PX 
QL = SORT (CX**2 + UY**2 + GZ#*2) 

W(l,l) = PX/PL 
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W ( 1 f 2 ) = PY/PL 
W ( 1 * 3) = PZ/PL 
W ( 2 » 1 } = OX/QL 
W ( 2 » 2 ) = QY/OL 
W ( 2 , 3 ) = CZ/DL 
W ( 3 » 1 ) = RX/RL 
W ( 3 » 2 ) = PY/RL 
V (3 ,3 ) = RZ/RL 
DO 10 J = 1 , 12 
CO 10 1 = 1,12 

io ;; (i,j) = o.o 

CALL EULEP ( E J( 1 , 1 ) , T, 3) 

CALL MULTB (W,T, 3,3,3, 3,3) 
DC 22 J=1 ,3 
Z(1,J) = T( 1 , J ) 

Z ( 5 , J 1 = 1(2, J) 

Z(9,J) = T (3, J ) 

JP3 = J+3 

Z ( 3,JP3 ) = T ( 1 , J ) 

Z ( 7, JP3 ) = T ( 3 , J ) 

22 Z(11,JP3) = T (2 , J ) 

CALL EULER ( E J ( 1 ,2 ) , T, 3 ) 

CALL MULTB (W,T, 3,3,3, 3,3) 
DO 24 J = 1 ,3 
JP6 = J +6 

Z( 2 , JP6 ) = T ( 1 , J ) 

Z( fc , JPb ) - T (2 , J ) 

Z(I0,JP6) = T (3, J ) 

JP9 = J+9 

Z( 4,JPP) = 7(1, J) 

Z( E,JP9i = T (3, J ) 

24 Z (12,JP9) * 1 (2, J) 

&"TURN 

999 CALL ZZBOMfc (6HDC0S1B ,N‘tRROR) 
END 
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SUBROUTINE DCOS2 ( C K ,E J , Z,KC J ,KE J ,KZ ) 

DIMENSION C J ( KCJ , 1 ) » CJ(KFJ,li 5 Z(KZ,l) 

DIMENSION W(3,3), T(3,3) 

C 

C SUBROUTINE 10 CALCULATE FINITE ELEMENT-.. 

C DIRECTION COSINE MATRIX 

C FPF A COMP IN FD MEMBF ANE-PEND ING TRIANGLE PLATE ELEMENT. 

C THE DIRECTION COSINE MATRIX RELATES LOCAL COORDINATE DISPLACEMENTS 
C TO GLOBAL COORDINATE DISPLACEMENTS. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE PLATE TO LIE IN AN X-Y PLANE 
C WITH JOINT I AT THE ORIGIN, JOINT 2 LIES ALONG THE POSITIVE 
C X AXIS, AND JOINT 3 IS IN THE POSITIVE Y DIRECTION. 

C ROW ORDER (LOCAL COORDINATE ORDER) OF DIRECTION COSINE MATRIX IS 
C (DX,DY,TZ) JOINT 1, THEN JOINT 2, 3, NEXT 

C ( DZ ,TX , TY ) JOINT 1, THEN JOINT 2, 3 

C WHERE DX,DY,DZ ARE TRANSLATIONS AND TX,TY,TZ ARE ROTATIONS. 

C COLUMN ORDER (GLOBAL COORDINATE ORDER) OF DIRECTION COSINE MATRIX IS 
C (U,V,W,P,C,F) JOIN! 1, THEN JOINT : , 3. 

C WHERF -U,V,W ARE TRANSLATIONS AND P,C,R ARE ROTATIONS. 

C EULER ANGLE CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

C CALLS FORMA SUBROUTINES EULER, MULTB,ZZBCMB. 

C DEVELOPED E v WA BENNIE LD . FEBRUARY 1973. 

C LAST REVISION BY WA 6E.NF IELD . MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS 

C CJ = INPUT MATRIX OF GLOBAL X,Y,Z COORDINATES AT TRIANGLE JOINTS. 

C ROWS 1,2,3 CORRESPOND TO X.Y,Z COORDINATES. 

C COLS 1,2,3 CORRESPOND TO JOINTS 1,2,3. SIZE(3,3). 

C EJ = INPUT MATRIX OF EULER ANGLES (DEGRESS) AT TRIANGLE JOINTS. 

C ROWS 1,2,2 CORRESPOND TO GLOBAL X,Y,Z PERMUTATION. 

C COLS 1,2,3 CO°°ESPOND TO JOINTS 1,2,3. SIZE<3,3). 


C 

z 

= OUTPUT 

DIKE 

CTICN COSINE 

MATRIX. SIZE ( 16,18) . 

c 

KCJ 

= INPUT 

ADD 

DINE NS ION 

* 

CJ IN CALLING PROGRAM. 

c 

KEJ 

= INPUT 

ROW 

DIMENSION 

OF 

EJ IN CALLING PPOGRAM. 

c 

KZ 

= INPUT 

ROW 

DIMENSION 

OF 

Z IN CALLING PROGRAM. MIN=18. 


C 

C NERRCR EXPLANATION 

C 1 = DIMENSION SIZE LESS THAN IS. 

C 

NERROR=l 

IF (KZ .LT. 18) GO TO 999 
PX = CJU ,Z )-CJ(l,l ) 

PY = CJ (2,2 )— C J (2,1) 

PZ = CJ(2 ,2)-CJ(3,l ) 

PL = SQPT(PX**2 + PY**2 + PZ**2) 

P X = PY*(CJ(?, 3 )-CJ (3,1 )) - PZ*(CJ(2,3)-rj(2,l>) 

RY = PZ*(CJ(1 ,2)-CJ (1 ,1 )) - PX*(CJ(3,3)-CJ<3,1)) 

RZ = PX*(CJ(2,3)* ^J(2,l )) - PY*(CJ(1,3)-CJ(1»1>) 

RL = SORT ( RX^*2 + RY**2 + RZ**2! 

OX - RY*P 2 - PZ*PY 
CY = RZ*PX - PX*P2 
QZ - RX*P Y - FY*PX 
CL = SORT (GX**? + QY**2 + OZ**2) 

W(l,l) = PX/PL 
W( 1 .2) - PY/PL 
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W(l,3) = PZ/PL 
WC2.1) = CX/QL 
W ( 2 , 2 ) = QY/QL 
W' 2*3) = CZ/OL 
W ( 3 * 1 ) = RX/RL 
WC3 *2 ) = PY/PL 
W ( 3 *3 ) = PZ/RL 
DC 10 J=! ,18 
DC 10 1=1,18 
10 ZtI,J) = 0.0 
DP 50 MW=1»3 

CALL EULER ( FJ 1 1 ,NW) ,T ,3 ) 
CALL MULTF (W,T, 3,3,3, 3,3) 
IZZ = 3*(NW-1) 

JZZ = 6 * ( NW— 1 ) 

DP 50 JW= 1,3 

JZ = JZZ+JW 

Z(1ZZ-*- 1 , JZ ) = T ( 1 , JW ) 

Z1IZZ+ 2 , JZ ) = T ( 2 , JW ) 
Z(1ZZ+IC,JZ) = T1 3, JW ) 

JZ = JZ+3 

Z<IZZ+ 3, JZ ) = T ( 3 , JW ) 

Z(IZZ + 11,JZ) = T ( 1 , JW ) 

50 Z ( IZZ + 12 * JZ ) = T( 2, JW ) 

RETURN 

C 

999 CALL ZZBCMfc C 6HDCCS2 ,NERROR) 
END 
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SUE ROUT INC DCOS3C (CJ ,E J, Z,KC J,KE J ,KZ ! 

DIMENSION C J ( KC J , 1 ! * FJ(KEJ,1!, Z(KZ,11 
DIMENSION W(2 *3 ! « T(3,31 
C 

C SUBROUTINE TP CALCULATE FINITE ELEMENT... 

C DIRECTION COSINE MATRIX 

C FOR A RELTANGULAF SHEAR PANEL ELEMENT. 

C THE DIRECTION COSINE MATRIX RELATES LOCAL COORDINATE DISPLACEMENTS 
C TO GLOBAL COORDINATE DISPLACEMENTS. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE PANEL TO LIE IN AN X-Y PLANE 
C WITH JOINT 1 AT THE X-Y ORIGIN, JOINT 2 LIES ALONG THE POSITIVE 
C X AXIS* JOINT 3 IS IN THE POSITIVE X,Y DIRECTION, AND JOINT 4 LIES 
C ALONG THE POSITIVE Y AXIS. 

C ROW ORDER (LOCAL COORDINATE ORDER! OF DIRECTION COSINE MATRIX IS 
C DX1,DX2,DX?,DX4, DY1 , DY2 , DY3, DY4 

C WHERE DX,DY ARE TRANSLATIONS. 

C COLUMN ORDER CC-LOFAL COORDINATE ORDER! OF DIRECTION COSINE MATRIX IS 
C (U,V,W) JOINT 1, THFN JOINT 2, 3, 4. 

C WHERE U,V,W APE TRANSLATIONS. 

C EULER ANGLE CONVENTION IS GlCEAL X,Y,Z PERMUTA1 ION. 

C CALLS FORMA SUBROUTINES EULER , MULT6 .ZZBCMB . 

C DEVELOPED EY RL WOHLEN. APRIL 1974. 

C LAST REVISION BY WA BENFIFLD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS 

C CJ = INPUT MATRIX OF GLOBAL X,Y,Z COORDINATES AT PANEL JOINTS. 

C ROWS 1,2,3 CORRESPOND TO X,Y,Z COORDINATES. 

C COLS 1,2 ,3,4 CORRESPOND TO JOINTS 1,2,3,4. SI2E«3,4l. 

C EJ = INPUT MATRIX OF EULER ANGLES (DEGREES! AT PANEL JOINTS. 

C FOWS 1,2,3 CORRESPOND TO GLOPAL X,Y,2 PERMUTATION. 

C COLS 1,2 ,2,4 CORRESPOND TO JOINTS 1, 2,3,4. SJZE(3,4l. 

C Z = OUTPUT DIRECTION COSINE MATRIX. SI2E(F,121. 

C KC J = INPUT ROW DIMENSION OF CJ IN CALLING PROGRAM. 

C KEJ - INPUT ROW DIMENSION OF EJ IN CALLING PROGRAM. 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN-S. 

C 

C NEPRCR EXPLANATION 

C 1 = DIMENSION SIZE TCC SMALL . 

C 

NERROR=l 

IF (KZ .IT. £ ) CC‘ TO 999 
PX = CJ(1 ,? i-C J( 1,1 ) 

PY = CJ(?,?!-CJ(2,1! 

PZ = f J(3 ,2)-CJ(3,l 1 

PL = SOFT (PX* 1 *? + PY**2 ♦ PZ**2) 

CX = C J ( 1 ,4 ) — C J ( 1,1! 

CY = C J (2 ,4 )-CJ( 2 ,11 
C2 - C J(3,4)-f J(3,l ) 

CL = SCRT(CX**2 + CY* *2 + CZ**2) 

W(l,l) = PX/PL 
W ( 1 , 2 ! = PY/PL 
W ( 1 , 3 1 - P//PL 
W ( 2 , ! ) = CX/CL 
W (2,2) = CY/CL 
W ( 2 ,3 ) = GZ/OL 
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no 10 j=i » 12 

DO 10 1=1, P 
10 ZU.J) = 0.0 
DO 50 IJNT=1,<* 

CALL EULER <EJ(1,IJNT> ,T,3) 
CALL HL'LTF (W,T, 2,3,3, 2,3) 

JZZ = 3*( I JNT-Ll 

DO 50 JW=1 ,3 

JZ = JZ2+JW 

ZdJNT ,JZ) = T ( 1 / JW ) 

5C Z ( IJNT+4, JI) = T(2, JW ) 

RETURN 

999 CALL ZZBOMB C6KDC0S3C ,NERRCR ) 
• END 
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EULER 


SUERCUTINE EULER (E.R.KR) 

DIMENSION E(1),R(KR,1) 

C 

C C LCUL/'TF FULFR ANGLE ROTATION TRANSFORMATION MATRIX* 

C FU> ER ANGLE CONVENTION IS GLOBAL X T Y,Z PERMUTATION. 

C DEVELOPED EY C*“ EODLEY • MARCH 1973* 

C 

SUBROUTINE ARGUMENTS 

r = INPUT VECTOR OF JOINT EULER ANGLES (DEGREES). 

LOCATIONS 1»2»3 CORRESPOND TO THE GLOBAL X,Y t Z 
PERMUTATION. SIZE (3). 

R = OUTPUT EULER POTATION TRANSFORMATION MATRIX. SIZE(3«3). 

Kft * INPUT RON DIMENSION CF R IN CALLING PROGRAM. 

D1CR = ATAN2C1., 1.1/45. 

C 

Cl = CPS(E(l)*DTOP) 

C2 = CCS(E(2) *DTCP I 
C3 = CCS( E(3)*DT0R) 

51 = SIN(Etl) *DTCR ) 

52 = SIN(E(2l*DTOP) 

53 = SIN(E(3i*DTCR) 

C 

R ( 1 * 1 ) - C2*C3 

M!,2J = -C2*S3 
R 1 1 >3) = S? 

R (2 f 1 ) = C1*S3 + S1*S2*C3 

R(2,2) = C1*C3 - S1*S2*S3 

P ( 2 »3 ) = — S1*C2 
R (3 » 1 ) = S1*S3 - C1*S2*C3 

R ( 3 *2 J = S1*C3 •*- C1*S2*S3 

R(3,3) = CH-C2 


C 


RETURN 

END 
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SUBROUTINE FINEL ( XY Z * JDCF ,EUL .NUTEL »N J, 

* NUTM, NUTK, NUTLT, NUTST,NUTP,V,LV,KV, 

* KR X,KR J, KR E ,NUTMX,NUTK X ,NUT1 ,NUT2 ,NUT3 I 
DIMENSION XYZ(KRX,1), JD0E(KRJ,1), FUL(KRE» II » VIII « LVI1J 
DIMENSION Wl ( 2*- , 24) , W2(24,24), W3I24.24I 

DATA KW/24/, IBLANK/6H /, 11/1/ 

DATA NIT, NOT/5 ,6/ 

C 

C SUBROUTINE TO CALCULATE (ON GPTICN) FINITE ELEMENT... 

C ASSEMBLED MASS MATRIX ( ON NUTMI, 

C ASSFM5LFO STIFFNESS MATRIX (ON NUTK I , 

C ELEMENT LOCAL LOAD TRANSFORMATION MATRICES, IVECS (ON NUTLT I , 

C ELEMENT GLOBAL LOAD TRANSFORMATION MATRICES, IVECS (ON NUTXX) , 

C ELEMENT STRESS TRANSFORMATION MATRICES, IVfcCS (ON NUTSTJ, 

C ELEMENT UNIT LOAD BUCKLING MATRICES, IVECS (ON NUTB) • 

C IVEC GIVES ELFMFNT DOF INTO GLOBAL DOF. EXAMPLES... 

C IVEC (6 ) =E34 PLACES ELFMFNT DCF 6 INTO GLOBAL DOF 834. 

C IVEC ( 3 } =G OMITS ELEMENT DOF 2 FROM GLOBAL DOF. THIS CONSTRAINS 

C ELEMENT DCF 3 TO ZERO MOTION. 

C DATA ARRANGEMENT ON NOTH, NUTK FOR THE ASSEMBLED MATRICES IS IN 
C SPARSE ( Yl FORMA SUBROUTINE FORMAT. 

C DATA ARRANGEMENT ON NUTLT, NUTKX, NUTST, NUTB FOP EACH FINITE 
C ELEMENT (WRITTEN IN SUBROUTINE AXIAL, BAR, ETC I IS 
C WRITE ( NUTW ) NAMEW,NE L,NP. ,NC,NAMEL ( IBLANK, 1 = 1 ,5) , 

C ((W(I,J) ,1 = 1 ,NR ) , J=1,NC ) , { IVECi 1 1 , I=1»NC I 

C NAMFW = NAKELT,NAMEKX,NAMEST, OR NAMEB. 

C NAMEL = AXIAL, PAP, ETC . 

C LAST RECORD (TO DENOTE TERMINATION! IS, 

C WRITE (NUTW) IELANK,( II , 1=1 ,30) 

C THE FOLLOWING UTILITY TAPES USE BASIC FORTRAN READ, WRITE. DO NOT 
C USE THESE TAPES IN SPARSE (Y) FORMA SUBROUTINES WHICH USE FORMA 
C SUBROUTINES YIN, YOUT (BECAUSE THEY USE BUFFER IN, BUFFER CUT). 

C NUTLT, NUTST, NUTMX , NUTKX, NUTB. 

C THE FOLLOWING UTILITY TAPES USE FORMA YIN, YOUT. 

C NL'TM, NUTK, NUT1, NUT2, NUT 3. 

C CALLS FORMA SUBROUTINES AXIAL ,BAR .FLUID ,GRAVTY,PAGEHD*QUAD , 

C RECTSP ,TR*'C-L ,YRVAO? ,ZZbCMB. 

C DEVELOPED FY WA BENNIE LD , CS BODLEY, PL WOHLEN. JANUARY 1973. 

C LAST REVISION C Y PL WOHLEN. may 1976. 

C 

C **♦=,**=**** **********»♦**.****»************♦ ***♦****$*********- -,;*****♦*%$** 
C INPUT DATA READ IN THIS SUBROUTINE FROM NUTEL. IF NUTEL = 5, DATA IS 
C READ FROM CARDS. 

C 50 NAKEL FORMAT (A6) 


c 

IF 

(NAMEL .EG. fci-tRETURN) 

RETURN 





c 

IF 

(NAMEL .EG. 6HAXIAL 

) 

CALL 

AXIAL 

(SEE 

iUPPT 

FOR 

INPUT) 

c 

IF 

(NAMEL .EG. 6HE A p 

i 

CALL 

BAP 

(SEE 

SUBRT 

FOR 

INPUT) 

c 

IF 

(NAME L .FC. 6HFLUID 

) 

CALL 

FLUID 

(SEE 

SUBRT 

FOR 

INPUT) 

c 

IF 

(Name l .fo. 6HGpavty) 

CALL 

GPAVTY 

(SEE 

SUBRT 

FOR 

INPUT) 

c 

IF 

(NAMEL .FO. 6 H QUAD 

) 

CALL 

QUAD 

(SEE 

SUBPT 

FUR 

INPUT ) 

c 

IF 

(NAMEL ,tC. 6HRECTSP) 

CALL 

REO'^SP 

(SEE 

SUBRT 

FOP 

INPUT) 

c 

IF 

(NAMFL .FG. 6HTRNGL 

) 

CALL 

TRNGL 

(SEE 

SUBRT 

FOR 

INPUT) 

c 

GO 

TO 5 0 









C 

C DEFINITION OF INPUT VA e I A£ LE S . 
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NAMEL = AXIAL* BAR, ETC AS SHOWN ABOVE. GIVES SUBROUTINE CALLED. 

EXPLANATION OF INPUT FORMATS. NUMBER INDICATES CARD COLUMNS USED. 

A = ANY KEYPUNCH SYMBOL. 

X = CARO COLUMNS SKIPPED. 

slut************************************************************************ 
SUBROUTINE ARGUMENTS CALL INPUT I 

XYZ = MATRIX OF JOINT GLOBAL X,Y,Z LOCATIONS. ROWS CORRESPOND 
TO JOINT NUMBERS. CGLUMNS 1,2,3 CORRESPOND TO THE JOINT 
X , Y,Z LOCATIONS RESPECTI VELY. SIZE(NJ,3). 

MAY BE EiUIVALENCED TO V(ll IN CALLING PROGRAM. 

JDOF = MATRIX 0*- JOINT GLOBAL DEGREES OF FREEDOM. ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 
TRANSLATION DOFS AND COLUMNS *,5,6 CORRESPOND TO THE JOINT 
ROTATION DOFS. SIZE(NJ,6). 

MAY BE EC'U I VALE NC ED TO LVC1I IN CALLING PROGRAM. 

EUL = MATRIX OF JOINT EULER ANGLES (DEGRFES I • ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE 
•jLOPAL X ,Y. Z PERMUTATION. SIZE(NJ,3I. MAY BE 
ECUIVALENCED TO VCKRX*(XYZ COL DIMI+11 IN CALLING PROGRAM. 

NUTEL = LOGICAL NUMBER OF TAPE CONTAINING ELEMENT INPUT DATA FOR 

THIS SUBROUTINE AND SUBROUTINES AXIAL, ETC GIVEN BY NAMEL. 

IF NUTEL = 5, DATA WILL BE READ FROM CARDS. 

NJ = NUMBER OF JOINTS OR ROWS IN MATRT OFS (XYZ1, (JDOF), (EUL*. 

NUTM = LOGICAL NUMBER OF UTILITY TAPE ON WHICH ASSEMBLED 
MASS MATRIX IS OUTPUT IN SPARSE NOTATION. 

NUTM MAY BE ZERO IF MASS MATRIX IS NOT FORMED. 

USES FORMA YIN, YOUT. 

MUTK = LOGICAL NUMFER OF UTILITY TAPE ON WHICH ASSEMBLED 
STIFFNESS MATRIX IS OUTPUT IN SPARSE NOTATION. 

NUTK MAY BE ZERO IF STIFFNESS MATRIX IS NOT FORMEO. 

USES FORMA YIN, YCUT. 

NUTLT = LOGICAL NUMEER OF UTILITY TAPE (ft WHICH ELEMENT LO L 
LOAD TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. 

NUTLT MAY BF ZERO IF LOAD TRANSFORMATIONS ARE NOT FORMED. 

USES FORTRAN READ, WRITE. 

NUTST * LOGICAL NUMRER OF UTILITY TAPE ON WHICH ELEMENT 

STRESS TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. 

NUTST MAY PE ZERO IF STRESS TRANSFORMATIONS ARE NOT FORMED. 

USES FORTRAN READ, WRITF. 

NUTB = LOGICAL NUMEER OF UTILITY TAPE ON WHICH ELEMENT UNIT LOAD 
BUCKLING MATRICES AND IVECS ARE OUTPUT. 

NUTB MAY BE ZERO IF BUCKLING MATRICES ARE NOT FORMED. 

USES FORTRAN READ, WRITE. 

V = VECTOR WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

KRX = ROW DIMENSION OF XYZ IN CALLING PROGP AM. 

KRJ - FOW DIMENSION OF JDOF IN CALLING PROGRAM. 

KRE = ROW DIMENSION OF EUL IN CALLING PROGRAM. 

NUTMX = LOGICAL NUMBER OF UTILITY TARE ON WHICH ELEMENT 
MASS MATRICFS AND IVECS ARE STORED,. 

NUTMX MAY PE ZERO IF MASS MATRIX IS NOT FORMED. 

USES FORTRAN READ, WRITE. 
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NUTKX = LCGICAL DUMBER CF UTILITY TAPE ON WHICH ELEMENT 

STIFFNESS MATRICES (SAME AS GLOBAL LOADS TRANSFORMATION 
MATRICES) AND IVECS APE STORED. 

NUTKX MAY bE ZERO IF STIFFNESS MATRIX IS NOT FORMED. 
USES FORTRAN RE AD , WRITE. 

NUT1 = LOGICAL NUMBER CF UTILITY T.»PE. USES FORMA YIN* YOUT. 

NUT2 = LOGICAL NUMBER C'F UTILITY TAPE. USES FORMA YIN, YOUT. 

NUT3 = LOGICAL NUMBER OF UTILITY TAPE. USES FORMA YIN, YOUT. 

NFRRPR EXPLANATION 
1 = NAMEL IMPROPERLY DEFINED. 


1001 FORMAT 

2001 FORMAT 

2002 FORMAT 

2003 FORMAT 

* 

* 

♦ 

* 

* 

200 * 


SUBROUTINE 

SUBROUTINE 


( Afc ) 

(//41X 35HJ0INT DATA USED IN 
(//35X 47H JOINT DATA USED IN 
( /16X . 8HDEGR EE S OF FREEDOM 

lex ?eHGLOBAL CARTESIAN COORDINATES 
12X 22HEULER ANGLES (DEGREES) 

/l^X 11HTRANSLATI0N 8X 8HR0TATI0N 
/ 2X5HJOINT 6X1HU 5X1HV 5X1HW C X1HP 
I1X1HX 11XIHY 11X1HZ 14X1 HX 10XIHY 

- 4X 3F1I.4) 


FINEL) 

FINEL (CONTINUED)! 


FORMAT (IX 15, 

3X 

616, 3X 3 FI 2 

IF 

(NUTMX .GT. 

0) 

REWIND 

NUTMX 

IF 

(NUTKX .C-T. 

0) 

REWIND 

NUTKX 

IF 

(NUTE .C-T. 

0) 

REWIND 

NUTE 

IF 

4 NUT L T .GT. 

G) 

PEW IND 

NUTLT 

IF 

(NUTST .GT. 

0) 

REW IND 

NUTST 


5X1H0 5X1HR 
10X1HZ /! 


DETERMINE SIZE OF FINAL MASS-STIFFNESS MATRIX FROM THE MAXIMUM DOF 
NUMBER IN JDOF . 

NDOF = JDOF (1,1) 

DO 35 1=1, NJ 
DO 35 J=1 ,6 

IF (JDOF ( 1, J) .GT. NDOF) NDCF=JDOF ( I, J) 

35 CONTINUE 

PRINT JOINT DOF, XYZ COORDINATES, EULER ANGLES. 

CALL PAGFHD 
WRITF (NOT, 2001) 

WRITE (NOT, 2003) 

NLINF = 0 

00 40 I J= 1 »NJ 

NLINF = NL1NE+1 

IF (NLINE .LE. 42) GO TO 40 

CALL PAGFHD 

WRITE (NOT, 2002) 

WRITE (NOT, 2003) 

NLINF = 1 

40 WRITE (NOT, 2004) IJ, (JDCF(1J,J), J=l,6), (XYZ(IJ,J), J~l,3), 

* (EUL(IJyJ), J=l,3) 

READ FINITE ELEMENT TYPE. 

50 READ ( NUT EL ,1001 ) NaMEL 

IF (NAMEL .EQ. 6HRE TURN ) GO TO 50C 
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IF 

(NAM EL 

. FC . 6HAXIAL 

) 

GC 

TO 

110 

IF 

(NAM EL 

.EQ. 6HEAR 

) 

GO 

TO 

140 

IF 

(NAMEL 

• EQ. 6FTRNGL 

) 

GO 

TO 

150 

IF 

(NAMEL 

.EQ. 6H FLUID 

) 

GO 

TO 

151 

IF 

(NAMFL 

•EQ • 6HQUAC 

» 

GO- 

TO 

160 

IF 

(NAMEL 

.EQ. 6HPECTSP) 

GO 

TO 

162 

IF 

(NAMEL 

.EQ. 6HGPAVTY) 

GO 

TO 

171 


NERR0R=1 

GO TO 999 

C BAR FINITE ELEMENT (AXIAL ONLY). 

110 CALL AXIAL ( XY2 » JDCF ,EUL »NUTEL *NJ» 

* NUTMX, NUTKX, NUTLT,NUTST, 

* V.,W2,W3,KRX,KRJ,KRE,KW) 

GO TO 50 

C BAR FINITE ELEMENT (COMBINED AXIALf TORSION, BENDING). 

140 CALL BAR (XY2 , JDCF ,EUL ,NUTEL ,NJ, 

* NUTMX, NUTKX, NUTB ,NUTLT,NUTST, 

* W1 , W2, W3 , KRX ,KRJ , KPE »KW) 

GO TO 5C 

C TRIANGULAR PLATE ELEMENT. 

150 CALL TRNGL ( XYZ , JDCF ,EUL »NUTEL, NJ, 

* NUTM X , NU TKX , NUTB ,NUTLT ,NUTST, 

* W1 »W2,W3 ,KRX »KRJ ,KRE*KW) 

GO TO 50 

C FLUID ELFMENT. 

151 CALL FLUID ( XYZ, JDCF ,EUL ,NUTcL,NJ, 

* NUTMX,NUTKX, NIJTLT,NUT$T, 

* W1,W2,W3,KRX,KRJ,KRE,KW) 

GO TO 5C 

C QUADRILATERAL PLATE ELEMENT. 

160 CALL QUAD ( XYZ , JDCF ,EUL ,NUTEL ,NJ , 

* NUTMX,NU TKX, NUTB ,NUTLT ,NUTST, 

* Wi,W2,W3,KRX,KRvi,KRE,KW) 

GO TO 5C 

C RECTANGULAR SHEAR PANEL. 

162 CALL RECTSP ( XYZ , JDCF ,EUL , NUTEL ,NJ , 

* NUTMX,NUTKX ,NUTLT ,NUTST, 

* W1,W2,W3 ,KRX,KPJ,KRE,KW) 

GO TO 50 

C GRAVITY FLBMRNT. 

171 CALL GRAVTY ( XYZ , JDCF ,EUL ,NUTEL,NJ , 

* NUTKX, 

* W1,W2,W3,KRX,KRJ,KRE,KW) 

GO TO 50 

TERMINATE FINITE ELEMENT DATA ON STORAGE DISKS. 

500 IF ( NUTMX .GT. 0) WRITE ( NUTMX ) I PLANK, (II, 1=1, 30) 

IF (NUTKX .GT. t) WRITE (NUTKX) IBL ANK, ( I 1 , 1= 1 ,30 ) 

IF (NUTP .GT. 0) WPITE (NUTB) IBL ANK, ( 1 1 , 1= 1 , 30 ) 

IF ( NUT L T .GT. C) WRITE ( NUTLT ) IBL ANK, ( 1 1 , 1 = 1 ,30 ) 

IF (NUTST .GT. 0) WRITE (NUTST) IBL ANK, ( I 1 , 1=1 ,30 ) 

SUM FINITE ELEMENT MATRICES. 

IF (NUTM.GT •(.-) CALL Y ZERO ( NU TM ,NDC F ,ND OF ) 

IF (NUTK.GT.O) CALL Y ZERO (NUTK,NDCF,NDOF ) 
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IF (NUTMX .GT. 0) CALL YRVAD2 (NUTMX,NUTM,ND0F,W1 ,KW, V,LV,KV* 

* NUT1,NUT?,NUT3) 

IF (NUTKX .GT. 0) CALL YRVAD2 (NUTKX,NUTK ,NDOF ,W1 ,KW, V,LV,KV» 

* NUT1 »NUT2»NUT3 ) 

RETURN 

C 

999 CALL ZZBOMB (6HFINEL .NERRCR) 

END 
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SUbROUTIN f FLUID (XYZ ,JDCF, FUL, NUTEL, NJ , 

* NUTMX, NUTKX, NUTLT, NUTST, 

* W, T, S, KX, KJ, KE, KW) 

DIMENSION XYZ (KX ,1) ,JDOF(KJ,l ), EUL(KE,1), W(KW,1), 

* T ( KW , 1 ) , S(KW,1) 

DIMENSION CJ(3,8) ,EJ(3,8) ,IV1 (24 ) , I VTET ( 12 ) , JM (4, 16 ) , VL 1 10 1 * 


* 

DV(12)»DI ST (12*12) 

i , TV (24) 


DATA 

NR W , NR ST/24 » 1/, IBLNK/ 6 H 

/, 11 ST /O/ 

DATA 

NIT, NOT/ 5,6 / 



DATA 

NAMEL / 6 H FLUID / 



DATA 

kcj / ? /,;kjm / 4 / 



DATA 

KDIST / 12 / , I FB AD 

/ 1 / 


DATA 

JM/1,2,3,^, 3, 6,4,2, 

2 , 6 ,4,5, 

3', 5,1,2, 

* 

1,3, 6 , 5, 1,6, 4, 5, 

?4, 5 f 

1,2, 4,5, 

* 

4, 7, 8,5, 5,2, 7, 6 , 

4 ,2 ,3 9 7 1 

1,3, 8 , 6 , 

* 

1*6,8, 5, 1,3, 4, 8 , 

1,2,3,67 

8*?, 7 , 6 / 


C 

C SUBROUTINE to CALCULATE (ON OPTION/ FINITE ELEMENT ... 

C MASS MATRICES AND jtVECS (ON NUTMX i , 

C STIFFNESS MATRICES (SAME AS GLODAL LOAD TRANSFORMATION MATRICES? 

C AND IVECS (ON NUTKX), 

C PRESSURE TRANSFORMATION MATRICES AND IVECS (ON NUTST), 

C FOR FLUID ELEMENTS. 

C ELEMENT SHAPE KAY EE TETRAHEDRON, PENTAHEDRON, OR HEXAHEDRON. 

C MASS, STIFFNFSS MATRICES ARE IN GLOBAL COORDINATE DIRECTIONS. 

C GLOBAL COORDINATE ORDER IS 

C (U,V,W) JOINT 1 , CHEN JOINT 2 ,3,4, ( 5,6, 7, 8 ) . 

C WHERE U,V,W ARE TRANSLATIONS. 

C IVEC GIVES ELEMENT DOF INTO GLOBAL DOF. EXAMPLES... 

C IVEC (6 ) =834 PLACES ELEMENT DOF 6 INTO GLOBAL DOF 834. 

C I VEC ( 3 ) =0 OMITS ELEMENT DOF 3 FROM GLOBAL DOF. THIS CONSTRAINS 

C ELFMENT DCF 3 TO ZFRO MOTION. 

C PRESSURE TRANSFORMATION MATRIC.FS RELATE CHANGE IN PRESSURE (DUE TO 
C COMPRESSIBILITY) TO DEFLECTIONS IN THE GLOBAL COORDINATE DIRECTIONS. 

C PRESSURE CHANGE WITHIN THE FLUID ELFMENT IS CONSTANT. STATIC PRESSURE 
C DUF TO GRAVITY AND FLUID HEIGHT IS NOT INCLUDED. 

C DATA ARRANGEMENT ON NUTMX, NUTKX, NUTST FOR EACH FINITE ELEMENT IS 
C (W=M,K,ST) 

C WRITE (NUTWX) NAMEW,NEL,NP,NC,NAMEL,(IBLNK,I=1,5) , 

C ( ( W( I , J ) , 1= 1 » NR ) , J= 1 ,NC ) , (IVEC(I)»I =1,NC ) 

C CALLS FORMA SUBROUTINES TE GF PM , VCR OS S , VD OT ,ZZBOMB. 

C DEVELOPED BY C S bODLEY. FEBRUARY 1974. 

C LAST REVISION BY WA bENFIELD. MARCH 1976. 

C 

C INPUT DATA READ IN THIS SUbROUTINE FROM NUTEL. IF NUTEL = NIT, DATA IS 
C READ FROM CAROS. 

C N AM EM , NAM f:K »N AMELT , NAME ST FORMAT (4(A6,4X) 

C RC,EKM F OF MAT (2(5X,E10)) 

C 20 NFL, J1 » J2 »J3» J4,J5* J 6 *J7» JE FORMAT (915) 

C IF <J1 .EC. 0) RETURN 

C GO TO 20 

C 

C DEFINITION OF INPUT VARIABLES. 

C NAMEM = TYPE OF MASS MATRIX WANTED. 

C = Ml, LUMPED MASS MATRIX. 
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NAKEK 


NAMELT 
NAME ST 

RO 

BKM 

NEL 

J1 
J 12 
J3 
J4 


J5 

J6 


J7 

J8 


M2, QUAS I—l RROT AT 10NAL CONSISTENT MASS MATRIX. 

M3, IRROTATIONAL MASS MATRIX. 

= fcH OR 6HNPMASS, NO MASS MATRIX CALCULATED. 

TYP f OF STIFFNESS MATRIX WANTED. 

= KI, LINEAR DISPLACEMENT ASSUMED. 

= 6H OR 6HNOSTI F , NO ST IFFNF5S MATRIX CALCULATED. 

IDENTIFICATION NAME FOR LOAD TRANSFORMATION MATRICES. (NOT YET). 
IDENTIFICATION NAME FOR PRESSURE TRANSFORMATION MATRICES. 

= 6 H OR 6HN0STRS, NO PRESSURE TRANSFORMATIONS CALCULATED. 

MASS DENSITY. 

BULK MODULUS. 

FINITF ELEMENT NUMPEF. FOR REFERENCE ONIY, NOT USED IN 
CALCULATIONS. WRITTEN ON NUTMX, ETC. 

JOINT NUMBER AT ELEMENT VERTEX 1. 

JOINT NUMBER AT ELEMENT VERTEX 2. 

JOINT NUMBER AT ELEMENT VERTEX 3. 

JOINT NUMBEP AT FLEMENT VERTEX 4. 

FOR A TFTPAHEDR ON . FACE 1,2,3 MUST BE NUMBERED CLOCKWISE AS 
VIEWED FROM OUTSIDE THE ELEMENT. 

JOINT NUMBER AT ELEMENT VERTEX 5. (USED FOR PENTAHEDRON AND 
HEXAHEDRON) . 

JOINT NUMBER AT ELEMENT VERTEX 6. (USED FOR PENTAHEDRON AND 
HEXAHEDRON) . 

FOR A PENTAHEDRON, FACE 1,2,3 MUST BE NUMBERED CLOCKWISE AS 
VIEWED FROM OUTSIDE THE ELEMENT. FACE A, 5,6 IS NUMBERED IN 
THE SAME ORDER AS FACE 1,2,3. A LINE JOINING JOINTS I AND A 
MUST FORM AN EDGE OF THE PENTAHEDRON. 

JOINT NUMBER AT ELEMENT VERTEJ 7. (USED FOR HEXAHEDRON). 

JOINT NUMBER AT ELEMENT VERTEX 8. (USED FOR HEXAHEDRON). 

FOR A HEXAHEDRON, FACE 1,2, 3, A MUST BE NUMBERED CLOCKWISE 
AS VIEWED FROM OUTSIDE THE ELEMENT. FACE 5, 6,7,8 IS NUMBERED 
IN THE SAME ORDER AS FACE 1,2,3,A. A LINE JOINING JOINTS 1 
AND 5 MUST FORM AN EDGE OF THE HEXAHEDRON. 


EXPLANATION OF INPUT FORMATS. NUMBER INDICATES CARD COLUMNS USED. 

I = INTFGER DATA, RIGHT ADJUSTED. 

E = DECIMAL POINT DATA, ANYWHERE IN FIELD. EXPONENT RIGHT ADJUSTED. 
X = CARD COLUMNS SKIPPED. 


SUBROUTINE ARGUMENTS (ALL INPUT) 

XYZ = MATRIX OF JOINT GLOBAL X,Y,Z LOCATIONS. ROW!. CORRESPOND 
TO JOINT NUMBERS. COLUMNS ^,2,3 CORRESPOND TO THE JOINT 
X , Y ,Z LOCATIONS RESPECTIVELY. SIZF(NJ,3). 

JDOF = MATRIX OF JOINT GLOBAL DECREES OF FREEDOM. ROWS CORRESPOND 
TO JHTNT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 
TRANSLATION DOFS aND COLUMNS 4,5,6 CORRFSPOND TO THE JOINT 
ROTATION DOFS. SIZE(NJ,6). 

EUL - MATRIX OF JOINT EULER ANGLES (DEGREES). ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE 
GLOBAL X , Y , 7 PERMUTATION, SIZE(NJ,3). 

NUTEL - LOGICAL NUMBER OF TAPE CONTAINING ELEMENT INPUT DATA FOR 
THIS SUBROUTINE. IF NUT FI. = NIT, DATA IS READ FROM CARDS. 

NJ =• NUME-FR OF JOINTS OR ROWS IN MATRICES ( XY2 ), (JDOF), (EUL). 

NUTMX = LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 
MASS MATRICES AND IVFCS ARE OUTPUT. 
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NUTMX MAY PE ZERO IF MASS MATRIX IS NOT FORMED. 

USES FORTRAN READ, WRITE. 

NUTKX - LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 

STIFFNESS MATRICES (SAMBAS GLOBAL LOADS TRANSFORMATION 
MATRICES) AND IVECS ARE OUTPUT. 

NUTKX MAY BE ZERO IF STIFFNESS MATRIX IS NOT FORMED. 

USES FORTRAN READ, WRITE. 

NUTLT = LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT LOAD 

TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. (NOT YET). 

NUT ST = LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 

PRESSURE TRANSFORMATION MATRICES AND IVECS AR r " OUTPUT. 

NUTST MAY BE ZERO IF PRESSURE TRANSFORMATIONS ARE NOT-FORMED. 
USES FORTRAN READ, WRITE. 

W = MATRIX WORK SPACE. MIN SIZE(24,24). 

T = MATRIX WORK SPACE. MIN SIZE (24,24). 

S = MATRIX WORK SPACE. MIN SIZE(24,24). 

KX = ROW DIMENSION OF XYZ IN CALLING PROGRAM* 

KJ = ROW DIMENSION OF JDOF IN CALLING PROGRAM. 

KE - ROW DIMENSION OF EUL IN CALLING PROGRAM. c 

KW = ROW DIMENSION OF W, T, AND S IN CALLING PROGRAM. MIN=24. 

NERROR EXPLANATION 

1 = INCORRFCT TETRAHEDRON GEOMETRY. 

2 - INPUT JOINT NUMBER EXCEEDS MAXIMUM ALLOWABLE NUMBER CF JOINTS. 

3 = NUTMX NOT POSITIVE. 

4 - NUTKX MOT PGSITUVE. 

5 = NUTST NOT POSITIVE. 

1001 FORMAT (5 (A6,4X) ' 

1002 FORMAT ( 3 ( 5X , E 10 .0 , ) 

1003 FORMAT (915) 

2001 FORMAT (//25X 38H INPUT DATA FOR FLUID (TETRA, PENTA, OR 

* 21H HEXAHEDRON) ELEMENTS) 

2002 FORMAT (//20X 38H INPUT DATA FDR FLUID (TETRA, PENTA, OR 

* 33H HEXAHEDRON) ELEMENTS (CONTINUED)) 

2003 FORMAT ( /12X7HMAS S = A6,13X7HSTIF = A6,6X13HL0AD TRANS = AS, 

* 3X15HSTRESS TRANS = A6, 3X 

* / 15X,4hRC * lr 10.3 , 13X7HBULKM = E10.3, 

* // 9X7HE LEMENT 6 HJOINT 3 GX7HJ0TNT 2 6X7HJ0INT 3 6X7HJ0INT 4 

* fc/7HJ DINT 5- GX7HJ0INT h 6X7HJDINT 1 6X7HJ0INT 8 

* / 9X6HNUMBER ) 

2004 FORMAT ( 3X ,9( 8X , 1 5 ) ) 

3G01 FORMAT (51H * * * * * UNCONVENTIONAL JOINT NUMBERING *♦***,/ 

* 915 ) 


IF (II ST .EC. 1) GO TO 3 

1 1 ST = 1 

DO 4 1=1,4 

II = 3*1 - 2 

DO 4 J=l,4 

J1 = 3* J - 2 

4 CALL UNITY (D 1 ST ( 1 1 , J 1 ) ,3 ,KDI ST) 
C 

3 NLINF = 0 
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CALL PAGE HD 
WRITE (NOT, 2001) 

READ (N'J TEL, 1001 ) NAMEM ,N AMEK ,NAMELT,NAMEST 
READ (NUTEL*1G02 ) RO , EKM 

WRITE (WOT, 2003) NAMEM, NAMEK,NAMELT,NAMEST, 

* RO, BKM 
IF (NAMEM .NE. 6HM3 ) GO TO 20 
DO 2 1=1,12 

2 D I ST ( I , I ) = 2c 

20 READ ( NUTEL , 1 003 ) NEL * Jl, 12* J3» J4» J5* J6 , J7» J8 

NERR0R=1 

IF (J1.LE.0 .AND c IFRAD.EQ.-l) GO TO 990 
IF (Jl .LE. 0) RETURN 
NLINE = NLINE + 1 
IF (NLINE .LE, 42) GO TO 30 
CALL PAGEHD 
WRITE (NOT, 2002) 

WRITE (NOT, 2003) NAME M,NAMEK,NAMELT,NAMEST, 

* RO, BKM 
NLINE = 0 

30 WRITE (N0T,2004) NEL, J1,J2, J3,J4,J5,J6, J7?J6 

NERRC'R=2 

IF (Jl.GT.NJ .OR. J2.GT.NJ .OP. J3.GT.NJ .OR. J4.GT.NJ) GO TO 999 
IF (J5.GT.NJ .OR. J6.GT.NJ .OR. J7.GT.NJ .OR, J8.GT.NJ) GO TO 999 

FORM FINITE ELEMENT COORDINATE LOCATIONS * EULER ANGLES, REVADD IVEC. 

LR = 10 
NJN = 8 

IF ( J7 ,NE. 0) \yO TO 38 
LR = 6 
NJN = fc 

IF ( J5 .NE. 0) GO TO 38 
LR = 1 
NJN = 4 
C 

38 NCOL = 3* NJN 
DO 5 1=1, NCOL 
DO 5 J = I » NCOL 
W(I,J) = 0. 

S(I,J) = 0. 

5 T ( I » J ? = 0. 

c 

DO ao 1-1,3 

CJ( 1 ,1 ) = XY2 ( Jl ,1 ) 

C J ( 1 , 2 ) = XYZ(J2,I) 

CJ(I,? ) = XY7 ( J 3 , 1 1 
C J ( I ,4) = XYZ (J4,I) 

E J ( I , 1 ) = EUL ( Jl , 1 ) 

E J ( 1 ,2 ) = EUL ( J2 , I ) 

E J ( 1 ,3 ) = Fl'L ( J? , I ) 

E J ( I ,4 ) = EUL ( J4 , I ) 

I VI (1 ) = JDOF ( J 1 , I ) 

I VI ( 1+3 ) - JDOF ( J2, I ) 
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IV1 (1+6) = JD0F(j3, I) 

4C I VI ( 1+9 ) = JDCF(J4,I) 

IF (LR .EL. 1) GO TO 50 

DO 42 1=1,3 
C J ( I ,5 ) = XYZ ( J5 » I ) 

C J (1,6) = XYZ ( J6 , I ) 

E J ( I ,5 ) = EUL ( J5 » I ) 

E J ( I ,6 ) = EUL ( J6 , 1 ) 

I VI (1 + 12) = JDOF ( J5 ,1 ) 

42 IVK1+-15) = JDC F ( J6 , 1 ) 

IF (LP .EC. 6) GO TO 50 

DC 44 1 = 1 ,3 

C J (1,7) = XYZ(J7,I) 

C J (1,9) = XYZCJP.I) 

E J ( I ,7 ) = EUL(J7,1) 

EJ(i,e) = pul ( je,i) 

IV1 (1+18) = JDCF ( J7 »I ) 

':4 I VI ( 1+2 1 ) = J:TJF(Jfc t I) 

50 DO 52 L=1,LR 
LA = L 

IF ( LR. EG .10) LA=L+6 
DO 52 1=1 ,4 
JNO = JM( 1 ,LA ) 

LI = 3*1 - 2 
IVTET (LI ) = 3~JNC - 2 
I VTET ( L 1 + 1 1 = 3*JN0 - 1 
53 IVTET { L 1+2) = 3*JNC 

CALL TEGEOM (CJ,JM(1,LA ),VL(L),DV, KCJ.IFEAD) 

IF (IFBAD.NE.O) GO TO 51 

WRITE (NOT, 3001) NEL, J1 , J2* J3 , J4,J5 ,Jfc, J7,J8 
IFBAD = -1 

51 SK = R0*VL(L>/16.0 

IF (NAMEM .EQ. 6HM3 ) SM=RO*VL ( L )/20 .0 
IF (LR .GT. 1) SM = SM/2 . 

CALL REVADD ( 1 . ,DV, L, IVTE T, T, 7 ♦ 12 » LR , NCOL , 1 ,KW ) 

52 CALL REVADD SM,DIST, IVTET, IVTET, W, 12 ,1 2 ,NCOL,NCOL,KDIST,KW) 

IF (NAMEM .NE . 6HM1 ) GO TO 220 
^ DO 210 1 = 1 ,NCOL 
SAVE = O.C 
DO 215 J=1,NC0L 
SAVE = SAVE + W(I,J) 

215 WCI,J) = 0.0 
210 W(I,I) = SAVE 

220 IF (LR .EG. 1 ) GO TO 60 
DO 55 1=2, LR 
VL ( 1 ) = VL(1 ) + VL( I) 

DO 55 J = 1 ,NCOL 
55 T( I , J) = 7(i t j) ♦ T (I , J ) 

VL( 1 ) = VL(l)/2. 
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DC 56 J=2,NC0L 
56 T(1,J) = Ttl,J)/2. 

C 

60 DC 61 J=1,NCCL 

61 TV ( J ) = T (1 » J ) 

0C 65 J=1,NJN 
J1 = 3*J - 2 

65 CALL EULER ( E J (1 * J) ,S ( J1 , J1 ) ,KW) 

IF ( NAME ST .EO. 6H .OP. NAMEST .EQ. 6hN0STRS) GO TO 90 

CALL PRESS (C J,T t NJN,NCCL,KCJ,KW| 

CALL MULTA (T,S,>'JN, NCOL, NCOL, KW,KW) 

90 CALL PTAPA ( W ,S , NCOL, NCOL ,KW,KW I 
CALL MULTA (TV,S,1,NCCL,NCCL,1,KW) 

BOV = EKM/VLtU 

DO 70 1=1, NCOL 

DC 70 J=I,NCOL 

S(I,J) = BCV*TV(I)*TV(J) 

70 S(J,II = S(I,J> 

C 

IF (NAMEM .EC. t>H .OR. NAM EM ^.EQ. 6HNCMASS) GO TC 110 

NERRCR=3 

IF (NUTMX .LE. 0) GO TC 999 

WRITE (NUTMX) N/ MEM ,NEL ,NCCL, NCOL, NAMEL ,( IBLNK , 1=1,5 1 , 

* <(K(I,J), 1=1, NCOL), J=l, NCOL), ( IV1C I ) ,1=1 ,NCOL) 

C 

110 IF (NAMEK .EQ. 6h .OR. NAKEK .EQ. 6HN0STIF) GO TO 120 

NERR0R=4 

IF (NUTKX .LE. 0) GO TC 9»9 

WRITE (NUTKX) NA'IEK ,N EL , NCOL , NCOL, NAM EL , (IBLNK, 1=1, 5), 

* ((S(I,J), 1=1, NCCL),J=1, NCOL), (IV1( I),I=1*NCCL) 

C 

120 IF (NAMEST .EC. 6H .OP. NAMEST .EQ. 6HN0STRS) GO TO 2C 

NERR0R=5 

IF (NUT ST .LE. 0) GO TC 999 
NJNP1 = NJN ♦ 1 

CALL MULT (T , S , T (NJN PI , 1 ) ,NJN, NCOL , NCOL, KW,KW ) 

CALL MULTA (T,W, NJN, NCOL, NCOL, KW,KW) 

NR ST = 2* NJN 

WRITE (NUT5T) NAMEST, NFL, NR ST »NCOL ,NAME L, VL ( 1 ) , (I ELNK,I=1,4-|, 

* J (T( I , J) , 1=1 ,NRST) * J=l , NCOL l , ( IV1 (I ), 1=1, NCOL) 

GO TC 20 

C 

999 CALL ZZBOME ( 6HFLUID ,NEP RGR ) 

END 



GRAVTY — 1/ 4 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 
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SUBROUTINE GRAVTY (XYZ, JDOF ,EUL , NUTEL ,NJ, 

* NUTKX , 

* W, T, S, KX, KJ, KE, KW) 

DIMENSION XYZIKX, 1 ) , JDOF( KJ ,1 ) , EUUKE* II * W(KW,1), 

* (KW,1), S(KW,1) 

DIMENSION CJ (3,4) ,EJ( 3,4) ,1 VI ( 12) ,1 VTRI (9) , JM ( 3,4 ) ,GV(3 ) ,EV(3 ) 
DATA IBLNK/6H / 

DATA NIT, NOT/ 5,6 / 

DATA NAMEL / 6HGRAVTY / 

DATA KCJM / 3 / 

DATA JM / 1,2,3, 1,3,4, 1,2,4, 4,2,3 / 

SUBROUTINE TO CALCULATE ION OPTION) FINITE ELEMENT ... 

STIFFNESS MATRICES (SAME AS GLOBAL LOAD TRANSFORMATION MATRICES) 
AND IVECS (ON NUTKX), 

FOR GRAVITY ELEMENTS. 

STIFFNESS MATRICES AR* 9- GL02I- 300R49-1T AIR 090-S8 
GLOBAL COORDINATE ORDER IS 

(U,V,W) JOINT 1, THEN JOINT 2,3, (4). 

WHERE U,V,W ARE TRANSLATIONS. 

IVEC GIVES ELFMENT DOF INTO GLOBAL DCF. EXAMPLES... 

IVEC (6) =834 PLACES ELEMENT DOF 6 INTO GLOBAL DOF 834. 

IVEC (3 ) =0 OMITS ELEMENT DOF 3 FROM GLOBAL DOF. THIS CONSTRAINS 
ELFMENT DOF 3 TO ZERO MOTION. 

DATA ARRANGEMENT ON NUTKX FOR EACH FINITE ELEMENT IS 
(W=K ) 

WRITE (NUTWX) NAMEW ,NEL,NR,NC,NANEL,( I8LNK»I=1,5) , 

((W(I,J),I=1,NR ) * J=1 *NC ), (IVEC(I),I-1,NC) 

CALLS FORMA SUBROUTINES KGRAV ,MULTA ,MULTE ,VCRCSS, ZZBOMB . 

OE VELPPEO BY C S BPOLEY. FEeRUARY 1974. 

LAST REVISION BY WA BENFIELD. MARCH 1976. 


********** ********** ******************** ************ ************** ******** 


INPUT DATA RF AD IN THIS SUBROUTINE FROM NUTEL. 
READ FROM CARDS. 

NAMEM ,N AMEK 
RC 

(GV(I), 1*1,3) 

20 NEL.Jl.J? ,J3,J* 

IF (J1 ,EQ. 0) RETURN 

GC TO 20 


IF NUTEL = NIT, DATA IS 

FORMAT (2 ( A6,4X ) 

FORMAT (5X,E10) 

FORMAT (3(5X» E10) ) 
FORMAT (515) 


DEFINITION OF INPUT VARIABLES. 

NAMEM = TYPE CF MASS MATRIX WANTED. 

= fcH PR 6HN0MASS, NO MASS MATRIX CALCULATED. 

NAMEK = TYPF OF STIFFNESS MATRIX WANTED. 

= Kl, LINEAR DISPLACEMENT ASSUMED. 

= 6 H OR 6HN0STIF, NO STIFFNESS MATRIX CALCULATED. 

RO = MASS DENSITY. 

GV = GRAVITY VECTOR. 

NEL = FINITE ELEMENT NUMBFR. FOP P E FERENCE ONLY, NOT USED IN 
CALCULATIONS. WRITTEN ON NUTKX. 

J1 = JOINT NUMBER AT ELEMENT VERTEX 1. 

J2 = JOINT NUMFER AT ELEMENT VERTEX 2. 

J3 = JOINT NUMBER AT ELEMENT VERTEX 3. 
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J4 = JOINT NUMREP AT ELEMENT VERTEX 4. (USED FOR QUADRILATERAL). 

THE ELEMENT MUST BE NUMBERED CLOCKWISE AS VIEWED FROM THE FLUID SIDE 
OF THE ELEMENT. 

EXPLANATION OF INPUT FORMATS. NUMBER INDICATES CARD COLUMNS USED. 

I = INTEGEP DATA, FIGHT ADJUSTED. 

F = DECIMAL POINT DATA, ANYWHERE IN FIELO. EXPONENT RIGHT ADJUSTED. 

X = CARD COLUMNS SKIPPED. 

********************** ****************************** ********************** 

SUBROUTINE ARGUMENTS CALL INPUT) 

XYZ = MATRIX OF JOINT GLOEAL X,Y,Z LOCATIONS. ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 
X ,Y ,Z LOCATIONS RESPECTIVELY. SIZE(NJ,3). 

JDOF = MATRIX OF JOINT GLOEAL DEGREES OF FREEDOM. ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 
TRANSLATION DOFS AND COLUMNS 4,5,6 CORRESPOND TO THE JOINT 
ROTATION DOFS. SIZE (NJ,6) • 

EUL = MAT FIX OF JOINT EULER ANGLES {DEGREES ) • ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE 
GLOEAL X,Y,2 PERMUTATION. SIZECNJ,3). 

NUTEL = LOGICAL NUMFE r OF TAPE CONTAINING ELEMENT INPUT DATA FOP 
THIS SUBROUTINE. IF NUTEL = NIT, DATA IS READ FROM CARDS. 

NJ = NUMEEP OF JOINTS OR ROWS IN MATRICES (XYZ), (JCOF), (EUL). 

NUTKX = LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 

STIFFNESS MATRICES (SpME AS GLOBAL LOADS TRANSFORMATION 
MATRICES) AND IVECS ARE OUTPUT. 

NUTKX MAY PE ZERO IF STIFFNESS MATRIX IS NOT FORMED. 

USES fortran read, writf. 

W = MATRIX WORK SPACE. MIN SIZE(12,12). 

T = MATRIX WORK SPACE. MIN SIZE(12,12). 

S = MATRIX WORK SPACE. MIN SIZE (12,12). 

KX = ROW DIMENSION OF XYZ IN CALLING PROGRAM. 

KJ = POW DIMENSION OF JDOF IN CALLING PPOGRAM. 

KE = ROW DIMENSION OF FUL IN CALLING PROGRAM. 

KW = ROW DIMENSION OF W, T, AND S IN CALLING PROGRAM. MIN=12. 

NFRRO° EXPLANATION 

1 = INPUT JOINT NUMBER EXCEEDS MAXIMUM ALLOWABLE NUMBER OF JOINTS. 

2 = NUTKX NOT POSITIVE. 

1001 FORMAT ( 5 (A6,4X ) ) 

100? FORMAT (3 (5X,E 10.0) ) 

1003 FORMAT (5!5) 

2001 FORMAT ( / /2 5 X ^ EH INPUT DATA FOR GRAVITY STIFFNESS (TRIANGLE OR 

* ?4H QUADRILATERAL) ELEMENTS) 

2002 FORMAT (//2CX 45HINPUT DATA FOR GRAVITY STIFFNESS (TRIANGLE OR 

* 36F QUADR 1LATFR AL » ELEMFNTS (CONTINUED)) 

2003 FORMAT ( / 1 2X7HMAS S = A6,13X7HSTIF = A6,oX 

* / 15X,4HRC< = E10.3, 13 X5HGVX = E1C.3, 13X5HGVY = E10.3, 

* 13X5HGVZ = E10.3, 

* //15X7HFLFMENT 13X7HJOINT 1 13X7HJ0INT 2 13X7HJ0INT 3 

* 13X7F JOINT 4 

* /15X6HNUMPER) 

20C4 FORMAT ( 1 BX,9 ( 15 * 15X) ) 
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NLINE = 0 

CALL PAGEHD 

WRITE I NOT ,2001 1 

READ (Nt»TEL,1001 > NAM EM »NAMEK 

READ CNUTEL *1 0C2 ) RO 

READ (NUTFL, 10C2 ) (GV( I 1 ,1=1 ,3) 

WRITE (NOT *2003 1 NAME M*NAMEK, 

* RC, (GVC 11*1-1*31 

20 READ (NUTEL, 10031 NEL «J1, J2,J3,J4 
IF (J1 .LE. 0) RETURN 
NLINE = NLINE ♦ 1 
IF (NLINE .LE. 42) GO TO 30 
CALL PAGEHD 
WRITE (NOT, 20021 
WRITE (NOT, 20031 NAME M, NAME*, 

* RO, (GV( 11*1=1,31 
NLINE = 0 

30 WRITE (NOT, 20041 NEL, JI , J2, J3 * J4 

NERR0R=1 

IF (Jl.GT.NJ .CR. J2.GT.NJ .OR. J3.GT.NJ .OP. J4.GT.NJ1 GO TO 999 

FORM FINITE ELEMENT COORDINATE LOCATIONS, EULER ANGLES, REVADD IVEC. 

LP = 4 
NJN = 4 

IF (JA .NE. 01 GO TO 3£ 

LP - 1 
NJN = 3 

38 NCOL = 3»N^N 
DO 5 1=1, NCOL 
DO 5 J=1 * NCOL 

W(I*J1 = c. 

S(I,J1 -- G. 

5 T( 1 , J 1 = 0. 

DO 40 1=1,3 
CJ(I,1 ) = XY2 ( Jl , I 1 
C J ( I , 2 ) = XYZ(J2,I) 

C J (1,3) = XYZ ( J3, I 1 
E J ( T , 1 1 ' ntL(Jlfl) 

E J ( I ,? ) = EUL ( J2 , I ) 

E J ( 1 ,3 ) = t L f L ( J 3 , 1 1 
I VI ( I I = JDOF ( JI , I 1 
I VI (1+3 ) = JDOF (J2, 11 
40 iviii+ 6 ) = jocF(j:->,n 

IF (LP .EL. 1 ) GO TO 50 

DO 42 2=1,3 
C J ( I ,4 ) r XYZ ( J4 , 1 ) 

EJ(I,4) = E L'L ( J4 t I ) 

42 IVKI+9) = jrOE(J4,I) 
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50 G a S0RT<GVm**2 ♦ GV!2)**2 ♦ GV I 3 1**2 ) 

DO 51 1=1*3 

51 EVCI) = -GVII1/G 

DO 52 L-ULR 

CALL KGRAV (CJ, JM( 1, LI »EV, A,W,KW,KCJM1 
C 

DO 53 1=1,3 
JNO = JM(ItL) 

LI = 3*1 - 2 
IVTRKL! 1 = 3* JNO - 2 
I VTRKL1 + 1 1 = 3* JNO - 1 
53 IVTRHL1 + 2) = 3*JNC 
SS = R0*G*A/2^. 

IF I LR .GT. 1) S$ = SS/2. 

52 CALL REVADD I SS*W*IVTRI *IVTRI *S*9*9 t NC0L*NC0L*KW»KW J 
C 

DP 65 J=!,NJN 
J1 = 3* J - 2 

65 CALL EULER C F J f 1 , J 1 , T( J1 , J1 1 ,KW) 

CALL ETABA < S ,T,NCCL ,NCOL*KW ,KW) 

C 

IF CNAMEK .EG. 6H .CR. NAMEK .EC. 6HN0STIF) GO TO 20 

NERRCR=2 

IF fNUTKX .LE. 01 GT TO 999 

WRITE ( N’UTKX 1 NAMEK ,NFL ,NCPL ,NCOL *NAMFL * CIBLNK* 1 = 1 ,5 1 * 

* MSU.Jlt I=1»NCQL )*J=1»NC0L)» ( IV1C If *1=1 *NCOL) 

C 

GO TP 20 
C 

999 CALL ZZBOMP (6KGRAVTY.N ERROR! 

END 
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SUBROUTINE K1A1 ( Al ,A2,RL ,E «Z,TS , KZ ,KTS ) 

DIMENSION Z(KZ,1 ) « TS|KTS,1) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C STIFFNESS MATRIX, 

C STRESS TRANSFORMATION MATRIX, 

C FOR AN AXIAL PCD ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C ROD MAY RE LINEARLY TAPERED OR UNIFORM. 

C CONSTANT FCRCF ASSUMED. 

C STIFFNESS MATRIX IS IN LOCAL COORDINATF SYSTEM. 

C STRESS TRANSFORMATION MATRIX P.ELATFS STRESS AT ROD ENDS IN LOCAL 
C COORDINATE SYSTEM TO DFFLFCT IONS IN LOCAL COORDINATE SYSTEM. 

C ROW ORDER IN STRESS TRANSFORMATION MATRIX IS 
C SIGMA-X1, SIGMA-X2 

C WHERE SIGMA IS NORMAL STRESS. 

C SX1C-), SX2 (>) IS TENSION. SX1(+), SX2(-| IS COMPRESSION. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE ROD TO LIE ALONG THE X AXIS 
C WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

C LOCAL COORDINATE ORDER IS 
C DX1,DX2 

C WHERE DX IS TRANSLATION. 

C DEVELOPED EY PL WOHLFN. SEPTEMBER 1972. 

C LAST REVISION BY RL WOHLEN- SEPTEMBER 1973. 

C 

C SUBROUTINE ARGUMENTS 

C Al = INPUT CROSS— SECT! ON AREA At ROD END 1. 

C A2 = INPUT CPOSS— StCTI ON AREA AT ROD END 2. 

C RL = INPUT ROD LENGTH. 

C E = INPUT YOUNGS MODULUS OF ELASTICITY. 

C Z OUTPUT STIFFNESS MATRIX. SIZE(2,2). 

C TS = OUTPUT STRESS TRANSFORMATION MATRIX. SIZE(2*2). 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=2. 

C KTS = INPUT COW DIMENSION OF TS IN CALLING PROGRAM. MIN=2. 

C 

S = A1*E/RL 
R = A2/A1 

IF (APS (R— 1.) .GT. .01) S = ( A2— Al ) *E / ( RL*ALCG( R ) ) 

C STIFFNESS MATRIX. 

Z(l,l ) = S 
Z ( 1 ,2 ) =-S 
Z ( 2 , 1 ) =-* 

Z ( 2 ,2 ) = S 

STRESS TRANSFORMATION MATRIX. 

TSf 1,1 ) = 7(1,1 )/Al 
TS ( 1 ,2 ) = Z ( 1 , ? ) / A 1 
TS ( 2 , 1 ) = 2(2,1) /A2 
TS ( 2 , 2 ) = Z(2,2)/A2 
C 

RETURN 

END 
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SUBROUTINE KIP] t BI 1 ,P 1 2, Cl ,C2,A1 ,A2, SF ,RL, E ,G, Z*T$,KZ»KTS) 
DIMENSION 2<K2 t n t TS(KTS,1) 

DATA EPS/1. E-15/ 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C STIFFNESS MATRIX, 

C STRESS TRANSFORMATION MATRIX, 

C FOR A BENDING (PLUS SHEAR] PEAM ELEMENT WITH UNRESTRAINED BOUNDARIES. 
C BEAM MAY BE LINEARLY TAPERED OR UNIFORM. 

C UNIFORM SHEAR AND LINEAR PENDING MOMENT VARIATION IS ASSUMED. 

C SHEAR STIFFNESS USES SF*Al*G AND SF*A2*C . IF ANY OF THESE VARIABLES 
C ARE ZERO, THERE IS NO SHEAR DEFORMATION IN BENDING. 

C STIFFNESS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C STRESS TRANSFORMATION MATRIX RELATES STRESS AT BEAM ENDS IN LOCAL 
C COORDINATE SYSTEM TO DEFLECTIONS IN LOCAL COORDINATE SYSTEM. 

C ROW ORDER IN STRESS TRANSFORMATION MATRIX IS 
C T AU—X] , TAU— X2 *S IGMA— X 1 , SIG?’ A— X2 

C WHERE SIGMA IS NORMAL STRESS (MC/Il AND TAU IS SHEAR STRESS (P/AI. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE PEAM TO LIE IN THE X-Z PLANE 
C WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

C LOCAL COOPDINATE ORDER IS 
C DZ1,DZ2,TY1,TY2 

C WHERE DZ IS TRANSLATION AND TV IS ROTATION. 

C DEVELOPED PY RL WCHLEN. FEBRUARY 1973. 

C LAST REVISION BY RL WCHLEN. APRIL 1976. 

C 

C SUBROUTINE ARGUMENTS 

C BI1 = INPUT CROSS-SECTION AREA MOMENT OF INERTIA AT BEAM END 1. 

C BI2 = INPUT SAME AS BI1 AT FEAM END ?. 

C Cl = INPUT DISTANCE FROM PENDING NEUTRAL AXIS TO OUTER FIBER 

C AT BEAM END 1 . 

C C2 = INPUT SAME AS Cl AT EEAM END 2. 

C A1 = INPUT CROSS-SECTION AREA AT BEAM END 1. CAN BF ZERO FOR NO 

C SHEAR DEFORMATION IN BENDING. SHEAR STRESS IN STRESS 

C TRANSFORMATION WILL BF SET TO ZERO. 

C A2 = INPUT SAME AS A1 AT PEAM END 2. 

C SF =• INPUT SHAPF FACTOR (K > FOR SHEAR IN KAG. 

C USE SF=0.0 FOR NO SHEAR DEFORMATION IN BENDING. 

C SF= 1.0 FOR A SOLID CIRCULAR CYLINDER. 

C SF-.5 FOR A THIN WALLED CIRCULAR CYLINDER. 

C PL = INPUT ROD LENGTH. 

C E = INPUT YOUNGS MODULUS OF ELASTICITY. 

C G = INPUT SHFAR MODULUS OF ELASTICITY. CAN BF ZERO FOR NO SHEAR 

C DEFORMATION IN (ENDING. 

C Z = OUTPUT STIFFNESS MATRIX. SIZE(4,4). 

C TS = OUTPUT STRESS TRANSFORMATION MATRIX. SIZE (4,4). 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=4. 

C KTS = INPUT ROW DIMENSION OF TS IN CALLING PROGRAM. MIN=4. 

C 

C BENDING FLEX 1BIL ITY. 

PPI = EI2/PI1 
REIM1 = PEI-1. 

IF (ABSIRBIMI ) .LT. . 

F PR = F»PI1*P> IM1 
RBI LN = A LOG ( R P I ) 


01 ) GO TO 15 
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FI! = U5 - 1./R6IM1 ♦ RE1LN/RBIM1**2) * <RL**3) / EBR 
F 12 = Cl. - RBILN/RBIM1) * IRL**2) / EBR 
F 22 = PL* RPILN / EER 
GC TO 20 

15 FI 1 = RL**3 / {3.*E*EI1) 

F12 = P L**2 / (2.*E*BI1) 

F22 = RL/ <E*BI1) 

C SHEAR FLEXIBILITY. 

20 IF (SF.LT.EPS .OR. Al.LT.EPS .OR. A2.LT.EPS .OR. G.LT.EPS IGO TO 30 
RA = A2/A1 

IF ( ABS ( R A— 1 • ) .LT. .01) GO TO 25 - 

FI 1 = F 1 1 + RL * ALOGCRA) / ( SF*G*( A2-A1 )) 

GC TO 30 

25 FI 1 = FI 1 •*- RL/*SF*Al*G) 

C 

C BENDING + SHFAR STIFFNESS MATRIX. 

30 D = F11*F?2 - FI 2**2 
Z(l,l) = F22/D 
ZC1.2) =-Z(l»l) 

Z(i,3) =— F12/D 
Z ( I »4-) = <-RL*F22 ♦ F 12 1/D 
Z (2,2) = Z(l,l) 

Z(2,3) =— Z (1*3) 

Z ( 2 » A) =-Z!l,4) 

2(3.3) = FI 1/D 
Z(3,4) = (RL+F12 - F11J/D 

Z(4,4) = (F22*RL**2 - 2.*RL*F12 ♦ F11I/D 
C SYMMFTR I ZE LOWER HALF. 

DC 40 J-l»4 
DO 40 I=J,4 
40 Z( I,J) = Z( J, I ) 

C 

C STRESS TRANSFORMATION MATRIX. 

DO 55 J = 1 »4 

TS(1,J) = C.O 
TS ( 2 ♦ J ) = C .0 

IF ( A 1 .GT. 0.0) TS ( 1 » J ) = Z i 1 « J )/Al 
IF ( A2 .GT. 0.0) TS (2 » J ) = Z(2,J)/A2 
TS(3,J) = Z(3,J)*C1/BI1 
55 TS(4,J) = Z<4, J)*C2/E 12 
C 

RETURN 

END 
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SUBROUTINE K1C1 C TJ 1 , TJ2 ,R1 ,R2 ,RL ,G ,Z *TS ,KZ ,KTS ) 

DIMENSION Z(KZ,I>, 7S(KTS,1) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C STIFFNESS MATS IX , 

C STRESS TRANSFORMATION MATRIX, 

C FOR A TORSION ROD ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C ROD MAY PE LINEARLY TAPERED OR UNIFORM. 

C CONSTANT TOR DUE ASSUMFD. 

C STIFFNESS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C STRESS TRANSFORMATION MATRIX RELATES STRESS AT ROD ENDS IN LOCAL 
C COORDINATE SYSTEM TO ROTATIONS IN LOCAL COORDINATE SYSTEM. 

C ROW ORDER IN STRESS TRANSFORMATION MATRIX IS 
C TAU-X1,TAU-X2 

C WHERE TAU IS SHEAR STRESS. 

C STRESS IS + OR - AS RIGHT HAND AXIS BETWEEN END POINTS 1 AND 2. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE ROD TO LIE ALONG THE X AXIS 
C WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

C LOCAL COORDINATE ORDER IS 
C Tx 1 »TX? 

C WHERE TX IS ROTATION. 

C DEVELOPED BY RL WCHLEN. FEBRUARY 1973. 

C LAST REVISION BY RL WOHLEN. SEPTEMBER 1973. 

C 

C SUBROUTINE ARGUMENTS 

C TJ1 = INPUT CROSS-SECTION SAINT VENANTS TORSION CONSTANT CJJ IN JG 
C AT ROD END 1. E.G., TJ1= .5*PI*R!**4 FOR A SOLID CIRCULAR 

C CYLINDER. T Jl=2 .*PI*T*R1**3 FOR A THIN WALLED CIRCULAR 

C CYLINDER. 

C TJ2 = INPUT CROSS-SECTION SAINT VENANTS TORSION CONSTANT (J) IN JG 

C AT POD END 2. 

C R1 = INPUT DISTANCE FROM TORSION AXIS TO CUTER FIBER AT ROD END 1. 

C R2 = INPUT DISTANCE FROM TORSION AXIS TO OUTER FIBER AT ROD END 2. 

C RL =■ INPUT ROD LENGTH. 

C G INPUT SHEAR MODULUS OF ELASTICITY. 

C Z = OUTPUT STIFFNESS MATRIX. SIZE(2,2). 

C TS = OUTPUT STRESS TRANSFORMATION MATRIX. SIZE(2,2>. 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=2. 

C KTS = INPUT ROW DIMENSION OF TS IN CALLING PROGRAM. MIN=2. 

C 

S = TJ1*G/RL 
P = TJ2/TJ1 

IF ( APS (R — 1 • ) .GT. .Oil S = ( TJ2-TJ1 )*G / ( RL*ALCG(R ) ) 

C STIFFNESS MATRIX. 

Z(lfl) = S 
2(1,2} =-f 

zc?,n =-s 

Z ( 2 ,2 ) = S 

C 

C STRESS TRANSFORMATION MATRIX. 

TS (1 ,1 I = Z (1 ,1 )*R1/TJ1 
TS ( 1 , 2 ) = Z C 1 ,2 KR1/T Jl 
TS. (2,1 ) = Z (2,1 l*R2/TJ2 
TS ( 2 ,2 ) = Z (2,2 )*R2/TJ2 
C 



RETURN 

END 
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SUBPCUT I Nf K2A1 ( X2 ,X3 , Y3 ,TH, E , ANU « Z , T, P, KZ ,KT ,KR ) 

DIMENSION Z(KZ,1),T(KT, 1) ,R(KR,1) 

DIMENSION XE(3)»ET(3) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C STIFFNESS MATRIX, 

C STRESS TRANSFORMATION MATRIX, 

C FOR A MEMBRANE TRIANGLE PLATE ELEMENT WIT:. UNRESTRAINED BOUNDARIES. 

C QUADRATIC DISPLACEMENT (LINEAR STRAIN) FIELD IS USED. 

C STIFFNESS. MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C STRESS TRANSFORMATION MATRIX RELATES STRESS AT TRIANGLE VERTICES 
C IN LOCAL COORDINATE SYSTEM TO DEFLECTIONS IN THE LOCAL SYSTEM. 

C ROW ORDER IN STRESS TRANSFORMATION MATRIX IS 
C ( SIGMA— X , SIGMA— Y ,TAU— XY ) JOINT 1, THEN JOINT 2, 3. 

~ C WHERE SIGMA IS NORMAL STRESS AND TAU IS SHEAR STRESS. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE PLATE TO LIE IN AN X-Y PLANE 
C WITH JOINT 1 AT THE X-Y ORIGIN, JOINT 2 LIES ALONG THE POSITIVE 
C X AXIS, AND JOINT 3 IS IN THE POSITIVE Y DIRECTION. 

C LOCAL COORDINATE ORDER IS 
C (DX,DY,TZ) JOINT I, THEN JOINT 2, 3. 

C WHERF DX,DY ARF TRANSLATION. AND TZ IS ROTATION. 

C CALLS FORMA SUPRCUTINES BTAEA AND MULTA. 

C DEVELOPED BY CS EODLFY , WA B ENFIELD. MARCH 1973. 

C LAST REVISION BY CS BCDLEY. SEPTEMBER 1973. 

C 

C SUBROUTINE ARGUMENTS 

C X2 = INPUT LOCAL X COORDINATE LOCATION OF JOINT 2. 

C X3 = INPUT LOCAL X COORDINATE LOCATION OF JOINT 3. 

C Y3 = INPUT LOCAL Y COORDINATE LOCATION OF JOINT 3. 

C TH = INPUT PLATE THICKNESS. 

C E = INPUT YOUNGS MODULUS OF ELASTICITY. 

C ANU = I N C UT POISSONS RATIO. (E/2G)-1. 

C Z « OUTPUT STIFFNESS MATRIX. SIZE (9,9). 

C T - OUTPUT STRESS TRANSFORMATION MATRIX. SIZE(9,9). 

C R = INPUT MATRIX WORK SPACE. SIZF(8,9). 

C KZ = INPUT ROW DIMENSION OF l IN CALLING PROGRAM. MIN=9. 

C KT = INPUT ROW DIMENSION OF T IN CALLING PROGRAM. MIN=9. 

C KR = INPUT ROW DIMENSION OF R IN CALLING PROGRAM. MIN-B. 

C 

DO b 1=1 ,° 

DO 5 J=l,9 
5 T ( I , J ) = 0.0 
DO 10 1=) ,9 
DO 10 J=1 ,9 
10 Z( I,J) = 0.0 

IF (TH .LF. 0.0) RETURN 
X22 = X2*X 2 

Y3? = Y3*Y3 

X2Y3 = X2*Y?. 

SE1 = X3/X 2 
G = E / ( 2 . + ?.*ANU) 

DD = F*TH/(1. - ANU** 2) 

DNU = DD* ANU 
DG = G*T H 
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DO 15 1=1 ,8 
DC 15 J=l,9 
15 R ( I » J ) = 0. 

C 

FOO = X2Y3/2. 

F 10 = X?Y3*(1. + SE1J/6. 

F01 = X2Y3/6. 

F20 = X 2Y3* ( 1 . ♦ SE1 + SE1**2)/12. 

F 1 1 = X2Y?*(1. + 2.*$El)/24. 

F02 = X2Y3/12. 

C 

Z ( 1 » 1 ) = DD*FCC/X22 
Z ( 1 ,3 ) = DD*F 0 1/X 22 
Z ( 1 *6 ) = DNU*F00/X2Y? 

Z ( 1 ,8 ) = DNU*F 10/X2Y3 
Z(2 ,2) = DG+F00/Y32 
Z ( 2 ,3 ) = OG*F 1 0/Y32 
Z(2,4] i = 2.*DG*F01/Y32 
Z ( 2 » 5 J - PG*F00/X2Y3 
Z ( 2 *7 ) = 2.*DG*FI0/X2Y? 

Z(2,8) = DG*F 02 /X 2Y3 
Z ( 3 ,3 ) = DD*F02/X22 + DG*F20/Y32 
Z(3,4) = 2.*DG*F11/Y32 
Z (3,5) = DG*F10/X2Y3 
Z ( 3 * 6 ) = DNU*F01/X2Y3 
Z ( 3 ,7 ) = 2.*DG*F20/X2Y3 
Z (3 ,8) = DNP*F11/X2Y3 + DG*F11/X2Y3 
Z(4,4) = 4.*uG< F02/Y32 
Z (**»5 ) = 2.*DG*F01/X2Y3 
Z ( 4 ,7 ) = 4.*DG*F11/X2Y3 
Z (4 , 8 ) = ?.*0G*F02/X2Y3 
Z (5 *5 ) = DG*FOO/X 22 
Z<5,7) = 2.*DG*F10/X22 
Z ( 5 » 8 ) = DG*F 01/X22 
Z ( 6 , 6 ) = DD+F00/Y32 
Z ( 6 » 8 ) = DD*F 10/Y32 
Z ( 7 ,7 ) = 4.*DG*F20/X2 2 
Z ( 7 » 8 ) = 2 »*DG*F1 1/X2 2 
Z{8,8) = DD*F 20/Y32 + DG*F02/X22 
DO 20 1=1,8 
DO 20 J= I ,8 
20 Z(J,I) = Z( I , J ) 

C 

R(1 ,1 ) = -1. 

R(l,‘.) = 1 . 

R(2,l) = SE1 - 1. 

P ( 2 » 3 ) = -Y3 
P(?,M = -'FI 
R ( 2 ,7 ) = 1. 

R(?,«) = Y3 
R ( 3 , ? ) = Y? 

R { 3 » 6 ) = -Y? 

P <4, 3) = Y?*( 1 . - S t 1 ) 

R(4,6) = Y?*5 E 1 

R ( 4 , 9 ) = -Y3 



K2A 1 


— 3/ 3 


C 

C 


R ( 5 , 2 ) = -1. 

R(5,3) = X2 
R ( 5 » 5 ) - 1. 

R ( 5 ,6 ) = ~X2 
R ( 6 ,2 5 = i>E 1 - 1. 

R16.5) = —5 El 
R ( 6 » 6 ) = X? 

R ( 6 , £ ) = 1. 

R(6,«) = —X3 
R ( 7 , 3 ) = — X2 
P (7,6 ) = X2 

P (8,3) = X2* ( £ El ~ 1. ) 

R 1 1 , o ) - —X 3 
R ( 8 *9) = X2 

CALL BTABA ( Z , R, 8 ,9 ,KZ, KR ) 

Dll = DO/TH 
D 12 - ANU*D11 
033 = G 
X E ( 1 ) = 0. 

XE ( 2 ) = 1. 

XF ( 3 J = SE1 
E7 ( 1 ) = 0. 

ET( 2 ) = 0. 

ET( 3 ) = 1 . 

DC 30 1=1,3 
K1 = 3*1 - 2 
K2 = K1 + 1 
K3 = K 1 +2 
T ( K1 » 1 ) = D11/X2 
T ( K 1 , ?• > = D11*ET(I)/X2 
T { K 1 , 6 ) = D12/Y3 
T(K1 ,8) = P12*XF ( I)/Y3 
T (K2 , 1 ) = D12/X2 
T ( K 2 , 3 ) = D12*ET(1)/X2 
T ( K2 ,6 ) = D11/Y3 
T ( K2 , 8 ) s D 1 1 X F ( I ) /Y 3 
T ( K3 , 2 ) = D33/Y3 
T ( K3 , 3 ) = D33*XE (D/Y3 
T ( K3 , a ) = 2.*D33*ET (I }/Y3 
T (K 3 , 5 ) = D33/X2 
T t K3 ,7 ) = ?,*D33*XE (I )/X2 
30 T ( K 3 , 8 ) = ’*ET (1) /X 2 

CALL MULT. (T ,F ,9,8, 9,KT ,KR ) 

RETURN 

END 
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MieRC'UTIfcfc K2E 1 (X2,X3,Y3,TH,E,ANU,7,TS,T,KZ,KTS,KT) 

DIMENSION Z(KZ,?),TS(KTS,1 ) ,T(KT,1) 

DIMENSION R( 1C, 1C), IV EC (1C ) ,CCEF(9) ,XE(3),ET(3) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C ST1FCNESS MATRIX, 

C STRESS TRANSFORMATION MATRIX, 

C FDR A FENDING TRIANGLE PLATE ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C CUBIC DISPLACEMENT (LINEAR CURVATURE) FIELD IS USED. 

C STIFFNESS MATRIX IS IN LOCAL C CORD INATE SYSTEM. 

C STRESS TRANSFORMATION (aTRIX RELATES STRESS AT JOINTS IN LOCAL 
C COORDINATE SYSTEM TO DEFLECTIONS IN THE LOCAL SYSTEM. 

C ROW O c DF R IN STRESS TRANSFORMATION MATRIX IS 

C 'SIGMA-X, SIGMA— Y,TAU-XY) FOR (Z=TH/2) AT JOINT 1, THEN JOINT 2,3, 

C ( SIGMA-X, SIGMA-Y,TAU-XY) FQP <Z=-TH/2) AT JOINT 1, THEN JOINT 2,3. 

C WHERE SIGMA IS NORMAL STRESS AND TAU IS SH^AR STRESS. 

C TUE LOCAL CCORuINATE SYSTEM ASSUMES THE PLATE TO LIE IN AN X-Y PLANE 
C WITH JOINT 1 AT THE X-Y lRIC-IN, JOINT 2 LIES ALONG THE POSITIVE 
C X AXIS, AND JOINT 2 IS IN THE POSITIVE Y DIRECTION. 

c local crr c dinate order is 

C (DZ,TX,TY) JOINT 1, THEN JOINT 2, 3. 

C WHERE D7 IS TRANSLATION AND TX ,TY ARE ROTATIONS . 

C CALLS FORMA SLF c OUTINES FT At A AND MULTA. 

C DEVELOPED FY CS ECDLEY. MARCH 1973. 

C LAST REVISION R. Y OS EODLEY. SEPTEMBER 1973. 

C 

C SUBROUTINE ARGUMENTS 

C X2 = INPUT LOCAL X CC’ORD INATr LOCATION OF JOINT 2. 

C X3 - INPUT LOCAL X COORDINATE LOCATION OF JOINT 3. 

C Y3 = INPUT LOCAL Y COORDINATE LOCATION OF JOINT 3. 

C TH = INPUT PLATE THICKNESS. 

C E = INPUT YOUNGS MODULUS OF ELASTICITY. 

C ANU - INPUT POISSONS RATIO. (£/2G)-l. 

C Z = OUTPUT STIFRNESS MATRIX. S1ZE(P,9). 

C TS = OUTPUT L r CAL ST D ESS TRANSFORMATION MATRIX. SIZE(16,9). 

C T = INPUT MATRIX WORK SR.CE. S 1 7 E ( 10 , 1 C ) • 

C KZ = INPUT FOW DIMENSION OF Z IN CALLING PROGRAM. MIN=9. 

C KTS = INPUT FCW DIMENSION OF TS IN CALLING PROGRAM. MIN'=ie. 

C KT - INPUT ROW DIMENSION OF T IN CALLING PROGRAM. MIN^IO. 

C 

DO 10 1=1,9 

DO 10 J=!,V 

10 Z ( I , J ) = 0.0 

DC 11 1=1,10 

DC 11 J=l,° 

11 TS(IffJ) = C.C 

IF (TH .LE. 0.0) ‘-tTURN 
C 

DC 12 1=1,10 
DO 12 J = 1 , 1 0 

12 T( I t j ) = r . c 

X2? = X?*X? 

X?A = x??*x?? 

Y?-2 = Y * Y3 

Y3A = Y:2*Y3? 



SE1 = X3/X2 

SEE = 5 f 1*SE1 

SE3 = SE2*5E1 

SE^ = SE3*$E1 

SEC1 = tl . ♦ SEl )/3. 

SEC2 = SEC1**2 

SEC3 = SEC1**3 

G = t/(2. ♦ 2.*ANU> 

DD = 

ONU - OD*AN(J 

DG = (G*TH**3)/12. 

al = dd/x: 

BE = DNli/ (X22*Y32 ) 

GA = rD/Y?A 

DE = A. *DG/(X22*Y32 ) 


Tll.ll 

= 

i . 

T(2,3) 

=: 

!. 

T ( 3 » 2 ) 

= 

1. 

T (A»l ) 

- 

1 . 

T(a,2) 

= 

1. 

T|A,A) 


1. 

T|A,7) 

= 

1 . 

T(E,3I 


1. 

T ( B » B ) 

- 

K 

TIE, 8) 


1. 

T (6 ,2 ) 

= 

1. 

T(6,a) 


2. 

T(6 ,71 


3 . 

T(?,n 

- 

- • 

T ( 7,2 ) 

- 

SF1 

T(?»3 ) 


1. 

T (7 1 a ) 

- 

S r 2 

T(7,5) 

= 

5. El 

T( 7 ,6) 


1 . 

T(7,7) 

- 

SEE 

*: (7,3) 


S E 2 

T ( 7 »° ) 

- 

St 1 

T ( 7 » 1 0 ) 

= 1 . 

T { 8 ,3 ) 

= 

1 . 

T ( 8 * 5 J 

= 

sn 

T ( 8 ,fc ) 

- 

2 • 

TIB, 8) 


c r 2 

Tie, 9) 

- 

2.*SE! 

T (ft,10) 

- 3 . 

T I c , 2 ) 


1 . 

T(9,M 

- 

?.*SE 1 

7(9, B) 

- 

1 . 

7(9,7) 


?.*SF 2 

T 1 ° , 8 ) 

- 

2. *5 FI 

T ( ° »° ) 

- 

1 . 

T (10,1 ) 

= 1 . 

T (1C,?) 

= S F C 1 

T ( 1 C , 3 

) 

= !./:•• 

T (10, A ) 

= SEC 2 
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TC10,5) - SEC1/3. 

T C 10 ,6 1 = 1./0. 

TC1C,7) = 5CC3 
TC10, 8) = SFC2/3. 

T ( 1C***) = SEC1/9. 

TC1C.101 = 1./27. 

DO 5 1=1,10 
DO 7 J=I,1C 

7 R ( I , J ) = 0. 

5 RCI.I1 = 1. 

DC 100 L=1 , 10 
JBIG = 1 

A 1 = AEBCTfL,! II 

DC 15 J=2 ,10 
A2 = ABMTIUJ)) 

IF I A2 .LT. All GO 70 15 
A 1 = A? 

JEIG = J 
15 CCNTI*UE 

IVECCL) = JBIG 
ALJPIG = TCI, JBIG! 

DC 17 J=1,1C 
T C L , J ) = T CL, J l/ALJFI G 
17 R(L,J» = P CL,J)/ALJBIG 
DC 25 1=1,10 
A1JEIG = TCI, JBIG! 

IF (’ .EC. L) GC TC 2 5 
DC 3C J= 1 , 1C 

T C 1 , J ) = T ( T , J ) - A 1JEIG=T ( L, J ) 
30 R I I , J ) = F C I , J ) - A1JBIG*FCL, J) 
25 CCNT IF-UF 
100 CONTINUE 

DC 40 1 = 1 , 1C 
IP = IVFC C I ) 

DC ^0 J=1 ,IC 
40 T C IR , J ) = R ( 1 , J 1 
DC 50 1=1,10 
DC 5‘- 4=1,10 
50 MI.J) = T ( I , J ) 

DC 2C I = 1, IC- 
BC 1,2) = Y3*F Cl,') 
c C 1 , 3 ) = 

PCI, 5) = 

F f 1,6) = -X2*“ Cl ,M 

R C I ,8 ) = Y?.* c C I ,R ) 

20 RCI,9) = -XP^U.F-) 

CCBFC1) = 1./?. 

CCEFC?) = >?/lf. 

CCFFC?) = — ( X 2 ♦ X ? ) / 1 £• • 

C Pf F C 4 ) = 1./3. 
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CCFFI5) = Y3/1P. 

CCEFC6) = f?.*X2 - X3)/18. 

CCEF(7) - 1./?. 

CCEFC8) = -Y3/°. 

COEF(9) = I2.*X3 - X2J/16. 

DC 80 1=1,10 
DO 60 J=1 ,9 

80 R f I , J ) = R f I » J ) ♦ P.II ,10)*CDEF(J) 
C 

DC 55 1=1,10 
DO 55 J=1 ,10 
55 T(I»JI =0. 

C 

FOC = X ?*Y3/2 . 

F 10 = X2*Y3*I1. + SF1I/6. 

F01 = X2*Y5/6. 

F20 = X2*Y3*C1. ♦ SF1 ♦ SE2)/12. 
F1I = X2*Y3*(1. + 2.*SEII/24. 

F02 = X2*Y3/12. 

C 

T( 4 ,4 | = 4.*AL*FC0 
Tf4,fc) = 4.*EE*F00 
T|4,7| = 12.*AL*F10 
TCi.,6) = 4.*AL*F01 
T(4,9) = ^.*FE*F!0 
Tf4,10) = 12.*FE*F01 
T ( 5 ,5 ) = DF*F 00 
T( 5 *fi ) = 2.*DE*F10 
T(5,9) = 2.*DE*FCl 
7(6,6) = A.*GA*FCC 
T ( 6 ,7 ) = 12.*EE*F!0 
TC6,e) = 4. *8 F *FC1 
T(6,9) = 4.*GA*F10 
HE, JO) = 12.*GA*F01 
T ( 7,7 ) = 36 .* AL*F20 
T ( 7*8 ) = 12.*AL*Fll 
T 17,9) = 1?.*PF*F2Q 

T (7,10) = 3fc.*bE*Fll 

T ( 8 ,8 ) = 4.*aL*F 02 4.#D r ’i'F2C 

T { 8 ,9 ) = 4.*FE*F11 ♦ 4.*DE*F11 
T ( 8 , 10 ) = 12 . = F F*F02 
T ( 9 ,9 ) = 4.*GA=»F20 + 4.*DE*F02 
T (9,10) = 12 . *GA*F 1 1 
T (10,10 ) = ?6.*GA*FC2 
C 

DC 6C 1 = 1 ,10 
DC 60 J=I,1C 
60 T I J , 1 ) = T ( 1 , J ) 

CALL ETAF'A IT «R • 10,9, KT «1C) 

DC e5 1 = 1,*' 

DC 85 J = 1 , c 
85 ZII.jI = T 1 1 , J ) 

C 

DC 73 1=1 ,4 
DO 73 J = 1 ,10 



73 TCI.J) = O.C 

Dll = -6.*OD/CCX2*TH)**2) 

D21 = ANU*011 

D22 = -6.*0D/C CY3*TH)**21 

D12 = ANl**T22 

D33 = -12.*DG/CCX2*Y3)*TH**21 


XE ( 1 } = 0 
XEC2) =1. 

• 

XEC3) =S E 1 

ETI1) = 0 

* 

ETC 2 1 = 0 

• 

ETC?) *1. 
DO 75 1=1 

»3 

K1 = 3*1 

- 2 

K2 = K1 ♦ 

1 

K3 = XI ♦ 

2 

T CK1 »A ) = 

2. *011 

T(K1,6) = 

2. *012 

T CK1 »7) = 

6.*D11*XE Cl ) 

TCKl,e) = 

2.*D11*ETCI ) 

T(K1,0| = 

2.*D12*XE Cl ' 

TCK1 »1C ) = 

6.*012*ETCI 1 

T(K2,A) = 

2. *021 

T(K2,6) = 

2. *022 

T(K2,7) = 

6.*021*XECI 1 

T(K2,8) = 

2.*021*ETCI ) 

TCK?,9| = 

?.*D22*XE Cl ) 

TCK2ilO)= 

t>.*L22*ETCI ) 

T C K3»5 1 = 

033 

TCK3.8) = 

2.*D33*XECI 1 

75 TCK3.9) = 

2.*D33*ETCI ) 


CALL MULT A C T,R,S, 10 ,9,KT, 1C ) 
DC? 77 1 = 1 ,9 
IP9 = !♦*> 

DD 77 J = 1 ,9 
TSCI.J) = T C 1 * J ) 

77 TS(IP° t J) = — TS ( 1 » J 1 


RETURN 

END 
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SUBROUTINE K3C1 (X3, Y3 ,TH,G,Z,T,KZ ,KT) 

DIMENSION 2 (K2 » 1 ) * T(KT,1) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C STIFFNESS MATRIX, 

C STRESS TRANSFORMATION MATRIX, 

C FOR A RECTANGULAR SHEAR PANEL ELEMENT WITH UNRESTRAINED BOUNDARIES. 
C LINEAR DISPLACEMENT (CONSTANT STRAIN) FIFLD IS USED. 

C STIFFNESS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C STRESS TRANSFORMATION MATRIX RELATES PANEL SHEAR STRESS (CONSTANT) 

C IN LOCAL COORDINATE SYSTEM TO DEFLECTIONS IN THE LOCAL SYSTEM. 

C THF LOCAL COORDINATE SYSTEM ASSUMFS THF PANEL TO LIE IN AN X-Y PLANE 
C WITH JOINT 1 AT THE X-Y ORIGIN, JOINT 2 LIES ALONG THE POSITIVE 
C X AXIS, JOINT 3 IS IN THE POSITIVE X,Y DIRECTION, AND JOINT 4 LIES 
C ALONG THE POSITIVE Y AXIS. 

C LOCAL COORDINATE ORDER IS 
C DX1,DX2,DX3,D*<*, DY1, 0Y2, E Y3, OY4 

C WmERE DX,DY ARE TRANSLATIONS. 

L DEVELOPED BY ?L WOHLEN. APRIL 1974. 

C 

C SUBROUTINE ARGUMENTS 

C X3 = INPUT LOCAL X COORDINATE LOCATION OF JOINT 3. 

C Y3 = INPUT LOCAL Y COORDINATE LOCATION OF JOINT 3. 

C TH = INPUT PANEL THICKNESS. 

C G = INPUT SHEAR MODULUS OF ELASTICITY. 

C Z = OUTPUT STIFFNTSS MATRIX. SIZE (P, 6). 

C T = OUTPUT 5T C ESS TRANSFORMATION MATRIX. SIZE 11, 8). 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=£ • 

C KT = INPUT ROW DIMENSION OF T IN CALLING PROGRAM. MIN=I. 

C 

C STIFFNFSS MATRIX. 

C = TH*G/^. 

A = C*X3/Y3 

B = C*Y3/X3 

Z(1 ,1) = A 
Z ( 1 ,2 ) = A 
2(1,3) =— A 
zn,4) =-a 
Z ( 1 ,E ) = C 
Z(I,6) =-C 
2(1,7) =— C 
Z ( I ,8 ) = C 
2(2,2) = A 
Z ( 2 »3 ) =-A 
2(2, M = —A 
2(2,3) = C 
Z ( ? ,6 ) =-C 
2(2,7) =-C 
Z ( 2 ,R ) = C 
2(3,3) = A 
2 (3, a) = A 
Z ( 3 » 3 ) =-C 
2(3,0 - C 
Z( 3,7) = C 
Z ( 3 , f ) =-C 
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Z ( *+ ) — A 

ZIA,5J ~—C 
Z (A, 6) = C 
Z(A,7) = C 
Z(A,B) =“C 
Z ( 5 *5 ) = F 
Z(B, t) = -P 
2 f 5»7 ) =-E 
Z!5,e» = B 
Z (6*6) = F 
Z ( 6 ,7 ) = B 
2(6,8) = — E 
Z ( 7 , 7 ) = B 
Z(7,E) -~ p 
ZIP, 8) = P 

C SYMMETRIZE LCWER HALF. 

DO 1C J = l,8 
DC 10 I=J,P 
10 ZCI,J) = Z(J,I> 

STRESS TRANSFORMATION MATRIX. 
DC 20 J = l,f 

20 T( 1 , J ) = 2.*Z(3,J)/(TH*X3) 

RETURN 
FND 
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SUBROUTINE KGRAV (CJ,JM , E V,A,W ,KW, KC J) 

DIMENSION C J C KC J ,1>, JM( 1), EV(I), WIKW,1) 

DIMENSION F(3,^) ,R12( 3) ,R13 (3 ) ,VN (3 ) ,Ft 3,31 
DATA KEF / 3 / 

C 

C SUBROUTINE TO DERIVE STIFFNESS MATRIX FOR A TRIANGULAR GRAVITY ELEMENT. 
C CALLS FORMA SUBROUTINES MULTA ,MULT6 ,VCROSS. 

C DEVELCPFD FY C S EODLEY. FEBRUARY 1974. 

C LAST REVISION EY C S BOO LEY. NOVEMBER 1974. 

C SUBROUTINE ARGUMENTS 

C CJ = INPUT MATRIX OF JOINT COORDINATES. SIZE (3,4). 

C JM = INPUT VECTOR OF JOINTS DEFINING A TRIANGLE. SIZE 131. 

C EV = INPUT VECTOR NORMALIZED GRAVITY. SIZE = 3. 

C A = OUTPUT AREA. 

C W = OUTPUT STIFFNESS MATRIX. 

C KW = INPUT ROW DIMENSION SIZE OF W IN CALLING PROGRAM. MIN=9. 

C KC J = INPUT ROW DIMENSION OF CJ IN CALLING PROGRAM. MIN=3. 

C 

J1 = JMtl ) 

J2 = JM ( 2 ) 

J ? = JM ( 3 ) 

DO 5 1=1,9 
DO 5 J«T f 3 
W( I . J) - 0. 

5 E ( J ,1 ) - 0. 

DC 7 1=1,3 

P 12 ( I ) = CJ ( I , J2 ) - CJUtJl) 

R 1 3 ( I ) = CJ ( I » J3 ) - CJCI.JI) 

DC e J=l,3 
e F(i,j» - i . 

7 F|I,I) = 2 • 

C 

CALL VCPPSS (R12,R13,VN,VAMAG,VBMAG,A,SINAB) 

DO 10 1=1,3 
10 VN ( I ) = VN'(1)/A 
ACUM = 0.0 
DO IE 1*1,3 
II = 3*1 - 3 

ACUM = ACUM + VNII|*EVI!) 

DO IS J= 1 , 2 
E ( I , 1 1 + J ) = VMJ) 

IS W<I1+J,II = V N ( J ) 

CALL MOLT B ( F ,E ,3 , ? , °,KE F , KB F ) 

CALL MULTA <W , E , °,3 ,9 ,K W, KE F ) 

A - A*ACUM*ACUM*ACUM 
C 

RETURN 

END 
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SUBROUTINE Ml A I ( Al ,A2 , RL ,RO,Z , KZ ) 

DIMENSION 2 (KZ » 1 ) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C LUMPED MASS MATRIX 

C FDR AN AXIAL POD FLEMENT WITH UNRESTRAINED BOUNDARIES. 

C ROD MAY BE L INFA p LY TAPERED OR UNIFORM. 

C MASS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE ROD TO LIE ALONG THE X AXIS 
C WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

C LOCAL COORDINATE ORDER IS 
C DX 1 ,DX2 

C WHERE DX IS TRANSLATION. 

C DEVELOPED EY KL WOHLEN. JANUARY 1973. 

C LAST REVISION EY RL WOHLEN. SEPTEMBER 1973. 

C 

C SUBROUTINE ARGUMENTS 


c 

AI 

= 

INPUT 

CROSS— SECT ION 

AREA 

AT 

ROD 

END I. 

c 

A2 

= 

INPUT 

CRC'SS-SECTICN 

AREA 

AT 

ROD 

END 2. 

c 

RL 

= 

INPUT 

ROD LFNGTH. 





c 

RO 

= 

INPUT 

MASS DENSITY. 





c 

Z 


OUTPUT 

mass matrix. 

SIZE(2,2) 

• 


c 

KZ 

- 

INPUT 

ROW DIMENSION 

OF Z 

IN 

CALLING PROGRAM 


C 

W1 = Al*RL*RC/6. 

W2 = A2*PL*RC/6. 

7(1,1) = 2.*WI ♦ WZ 
Z ( 1 ,2 ) = 0. 

2(2,1) = e. 

Z (2 ,2 ) = rfl ■*- 2.*W2 
C 

RETURN 

END 
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SUE'S OUT I NO MIA? C A 1 ,A2»RL»R0,Z»KZ) 

DIMENSION ZIKZ,1) 

C 

C SUBROUTINE TC CALCULATE FINITE ELEMENT... 

C CONSISTENT MASS MATRIX 

C FOR AN AXIAL F CD ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C ROD MAY PE LINEARLY TAPERED CP UNIFORM. 

C LINEAR DISPLACEMENT FUNCTION ASSUMED. 

C MASS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE ROD TO LIE ALONG THE X AXIS 
C WITH JOINT 1 AT THE ORIGIN* JOINT 2 ALONG THE POSITIVE X AXIS. 

C LOCAL COORDINATE ORDER IS 
C DX 1 *DX2 

C WHFRE DX IS TRANSLATION. 

C DEVELOPED EY RL WCHLFN. SEPTEMBER 1972. 

C LAST REVISION BY RL WOHL EN . SEPTEMBER 1973. 

C 

C SUBROUTINE ARGUMENTS 


c 

Al 


INPUT 

CROSS— SECTION 

AREA 

AT 

ROD END 

1 . 

c 

A2 

= 

INPUT 

CROSS-SECTION 

AREA 

AT 

ROD END 

2. 

c 

RL 

=■ 

INPUT 

ROD LFNOTH . 





c 

RO 

= 

INPUT 

MASS DENSITY. 





c 

Z 

=: 

OUTPUT 

MASS MATRIX. 

SIZE ( 

2,2) 

• 


c 

KZ 

- 

INPUT 

ROW DIMENSION 

OF Z 

IN 

CALLING 

PROGRAM. 


C 


W1 

= A] 

i 

L*RC 

/I 

i* m 


W2 

= M 

'*R 

L* C C 

/I 

2. 


Z ( 1 

,1) 

= 

3 • 

1 

♦ 

W2 

Z C 1 

*2) 

- 

W 

1 

+ 

W? 

Z (2 

*1) 

- 

zu. 

2) 



Z ( 2 

*2) 

~ 

w 

1 

+ 

3.*W2 

RETURN 






END 
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SUBROUTINE M1P1 ( Al ,A2, RL,RO,Z,KZ ) 

DIMENSION Z(KZ,1) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C LUMPED MASS MATRIX 

C FOR A BENDING BEAM ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C BEAM MAY FE LINEARLY TAPERED OR UNIFORM. 

C MASS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE BEAM TO LIE IN THE X-Z PLANE 
C WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

C LOCAL COORDINATE ORDER IS 
C DZ1»DZ2»TY1 ,TY2 

C WHERE D2 IS TRANSLATION AND TY IS ROTATION. 

C DEVELOPED BY RL WCHLEN. FEBRUARY 1973. 

C LAST REVISION BY RL WCHL'N. SEPTEMBER 1973. 

C 

C SUBROUTINE ARGUMENTS 

C Al = INPUT CROSS— SFCT ION AR FA AT BEAM END 1. 

C A2 = INPUT CROSS— SECTION AREA AT BEAM END 2. 

C RL = INPUT E FAN! LENGTH. 

C RO = INPUT MASS DENSITY. 

C Z = OUTPUT MASS MATRIX. SIZE(4,4). 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=4. 

C 

W1 = Al*PL*RO/t . 

W2 = A2*RL*RO/fc. 

DO 10 J=1 ,4 
DO 10 1=1,4 


10 Z(I,J) 

= C.O 

Z (1 ,1 1 

= 2.*W1 + W 2 

Z (2,2) 

= W1 + 2.*W2 

Z ( 3 , 3 ) 

= ( Al*E.C*RL**3 )/24 

2(4,4) 

= (A2*R0*RL**3 )/24 

RETURN 


END 
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SUBROUTINE M1B2 ( A1 , A2» RL ,R0, 7 ,KZ ) 

DIMENSION Z(KZ,1) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C CONSISTENT MASS MATRIX 

C FOR A BENDING BEAM ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C BEAM MAY BE LINEARLY TAPERED OR UNIFORM. 

C CUBIC DISPLACEMENT FUNCTION ASSUMED. 

C MASS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE BEAM TO LIE IN THE X-Z PLANE 
C WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

C LOCAL COORDINATE ORDER IS 
C DZ1,DZ2,TY1 ,TY2 

C WHERE DZ IS TRANSLATION AND 7Y IS ROTATION. 

C DEVELOPED BY PL WOHLEN. FEBRUARY 1973. 

C LAST REVISION BY RL WOHLEN. SEPTEMBER '•973. 

C 

C SUBROUTINE ARGUMENTS 

C A1 = INPUT CROSS-SECTION AREA AT BEAM END 1. 

C A2 = INPUT CROSS-SECTION AREA AT BEAM END 2. 

C RL = INPUT BEAM LENGTH. 

C RO = INPUT MASS DENSITY. 

C Z = OUTPUT MASS MATRIX. SIZE(4,4). 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=4. 

C 

W1 = Al*PL*RC/E40. 

W2 = A2*RL*RC/840. 

RL? = RL**? 

Z(l,l) = 24C.*W1 + 72.*W2 

Z ( 1 ,2 ) = 54.*W1 + 5^.*K2 

Z ( 2 , 2 ) = 72 .*W1 + 240. *W2 

Z(l,3) 30.*W1 + 14.*W?)*PL 

Z(l,4) = ( 14.*W1 + 12.*W?)*RL 

Z ( 2 » 3 ) =-< 12.*W1 + 14.*W2)*RL 

Z ( 2 ,4) = ( 14.*W1 30.*W2)*RL 

Z (3 ,3 ) = ( 5.*W1 + 3 .*W2 )*PL2 

Z ( 3 ,4) =-( 3.*W1 + 3.*W2)*RL2 

Z(4,4) = ( 3.*wi ♦ 5.*W2)*RL2 

DO 10 J=1 ,4 
DO 10 1 = J ,4 
10 Z(I,J) = Z(J,I) 

C 

RETURN 

END 
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-DUB ROUTI N c MIC 1 ( PI 1, PI 2, RL ,RO,Z ,KZ ) 

DIMENSION Z (KZ » 1 ) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C LUMPED MASS MATRIX 

C FOR A TORSION ROD ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C ROD MAY PE LINEARLY TAPERED t'R UNIFORM. 

C MASS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE ROD TO LIE ALONG THE X AXIS 

C WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

C LOCAL COORDINATE ORDER IS 
C TX 1 , TX2 

C WHERE TX IS ROTATION. 

C DEVELOPED BY RL WOHLFN. FEBRUARY 1973. 

C LAST REVISION BY PL WCHLEN. SEPTEMBER 1973. 

C 

C SUBROUTINE ARGUMENTS 

C PI1 INPUT CROSS— SECTI ON POLAR AREA MOMENT OF INERTIA AT ROD END 1. 

C P 1 2 = INPUT CROSS-SECTION POLAR AREA MOMENT OF INERTIA AT ROD END 2. 

C RL = INPUT FDD LENGTH. 

C RO = INPUT MASS DENSITY. 

C Z OUTPUT MASS MATRIX. SIZE(2,2). 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=2. 

C 

W1 - PIl*PL*R0/6. 

W2 = ei?*RL*RC/6. 

ZO,l) = + W 2 

Z C 1 ,2 ) = 0. 

z(2,n = o. 

Z ( 2 , 2 ) = W1 + 2.*W2 
C 

RETURN 

END 



MIC? 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 


SUBROUTINE MIC,' (P1!,PI2,RL,R0,Z,KZ) 

DIMENSION Z(KZ t l) 

SUBROUTINE TO CALCULATE FINITE ELEMENT,.^ 

CONSISTENT MASS MATRIX 

FOR A TOP SION ROD ELEMENT WITH UNRESTRAINED BOUNDARIES. 

ROD MAY E'F LINEARLY TAPERED OR UNIFORM. 

LINEAR DISPLACEMENT FUNCTION ASSUMED. 

MASS MATRIX 15 IN LOCAL COORDINATE SYSTEM. 

THE LOCAL COORDINATE SYSTEM ASSUMES THE ROD TO LIE ALONG THE X AXIS 
WITH JOINT 1 AT THE ORIGIN, JOINT 2 ALONG THE POSITIVE X AXIS. 

LOCAL COORDINATE ORDER IS 
TX 1 »TX2 

WHERE TX IS ROTATION. 

DEVELOPED BY RL WCHLEN. FEBRUARY 1973. 

LAST REVISION BY PL WCHLEN. SEPTEMBER 1973. 


PI1 

PI2 

RL 

RO 

KZ 


UPROOT INE ARGUMENTS 

- INPUT CROSS-SECTION POLAR AREA MOMENT OF INERTIA aT ROD END 1. 

- INPUT CROSS-SECTION POLAR AREA MOMENT OF INERTIA AT ROD END 2. 
= INPUT POD LENGTH. 

= INPUT MASS DENSITY. 

= OUTPUT MASS MATRIX. SIZE (2,2). 

= INPUT ROW DIMENSION Of Z IN CALLING PROGRAM., MIN=2. 


W1 - PI!* 1 ?' *P 0/12 • 

W2 = PI2* f.L*RC>/12. 
Z(l,l) = 3.*Wl + W 2 

Z ( 1 , 2 ) = Vvl + W2 

Z (2 , 1 ) - 2(1,2) 

Z ( 2 , 2 ) = Wl + 3.*W2 

R FT URN 
END 
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SUP C CL'T INF M?Al (X2>V3,THrRC:*Z,K2) 

DIMENSION 2 ( K2 » 1 ) " \ 

C 

C SUBROUTINE T i' CALCULATE FINITE ELEMENT... 

C LUMPED MASS M h ’ IX, 

C FOR A MEMBRANE TRIANGLE PLATE ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C MASS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE PLATE TC LIE IN AN X-Y PLANE 
C WITH JOINT 1 AT THE X-Y ORIGIN, JOINT 2 LIES ALONG THE POSITIVE 
C X AXIS, AND JOINT 3 IS IN THE POSITIVE Y DIRECTION. 

C LOCAL COORDINATE OF OFF IS 
C CCX,DY,TZ> JOINT 1, THEN JOINT ?, 3. 

C .WHERE DXtDY ARE TRANSLATIONS AND TZ IS XTAT1CN. 

C DEVELOPED FY WA BBNF1ELD. FEBRUARY 1973. 

C LAST REVISION EY FL KGHLEN. SEPTEMBER 1973. 

C 

C SUBRCUT INt ARGUMENTS 

C X? = INPUT LOCAL > COORDINATE LOCATION CF JOINT 2. 

C Y3 - INPUT LOCAL Y COORDINATE LCCATICN CF JOINT 3. 

C TH = INPL ,T PLATE THICKNESS. 

C PC - INPUT MASS DENSITY. 

C Z = OUTPUT MASS MATRIX- SI2EC9,°). 

C KZ = INPUT ROW DIMENSION CF Z IN CALLING PROGRAM. MIN=9. 

C 

AP FA - C.5*X2*Y3 
CM = CR0*TF*ARrA)/3.0 
DC 1C 1 = 1, c 
DO 10 J = i, c 
10 ZJI,J) = C.C 
DO 2C 1=1, ° 

20 2(1,1) = CM 
RETURN 
END 



M2A2 


1/ 3 


SUP P PUT INI M2 A? <X2,X3,Y3,TH,PHC,Z,1,R*kZ,KT,IS.R.l 
DIMENSION UM,1), TOO, 1),F tKR,!) 

c 

C SUBROUTINE TP CALCULATE FINITE ELEMENT...' 

C C°\'SIST ,r NT MASS M ATR 1 1 X, 

C FOR A MEMBRANE TRIANGLE PLaTE ELEMENT WITH L’NRE STRAINED BOUNDARIES. 

C QUADRATIC DISPLACEMENT (LINEAR STRAIN) FIELD IS USED. 

C MASS MATRIX 15 IN LOCAL CC'CRDINaTE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM A c "‘MLS THF PLATE TO LIE IN AN X-Y PLANE 
C WITH JC1NT 1 AT THE X-Y CFIGIw, JOINT 2 LIES ALONG THE POSITIVE 
L X AXIS, AMD Jtlk; 3 IS IN THE POSITIVE Y DIRECTION. 

C LOCAL COORDINATE ORDER IS 
C inv,r.Y,TZ) JOINT 1, THEN JOINT 2, 3. 

C WHERE DX,DY A c E TRANSLATIONS AND TZ IS ROTATION. 

C CALLS PCPMA SUBROUTINES ETAtA. 

C DEVELOPED FY CS EOPLEY. M.ARCh 1®?3. 

C LAST REVISION FY HA EENF1ELD. SEPTEMBER 1°7--. 

C 

C SUBROUTINE ARGUMENTS. 

C X2 = INPUT LOCAL X COORDINATE LOCATION OF JOINT 2. 

C X3 -= INPUT LOCAL X CCCPDINATE LOCATION OF JOINT 3. 

C Y3 = INPUT LOCAL Y COORDINATE LOCATION OF JOINT 3. 

C TH = INPUT PLATE THICKNESS. 

C RHC - INPUT MASS TENSITY. 

C Z = OUTPUT MASS MATRIX. SI2E(9,9). 

C T = INP”T MATRIX WORK IPaCF. SI2F(10,10I. 

CP - INPUT MATRIX R'CF K. SPaCE. SIZE! 10,10). 

C KZ = INPUT FOW DIMENSION OF Z IN CALLING PROGRAM. KIN=9. 

C KT = INPUT F.CW DIMENSION OF T IN CALLING PPCGRAM. MIN=10. 

C KR = INPUT ROW DIMENSION OF P ] I_N CALLING PROGRAM. NIN=10. 

c - 

SE? = X3/X2 
SE? = Stl*SEl 
SE? - SENSEI 
SE** = SFi*SE 1 
X2Y3 = X2*Y2*RHO*TK 
C 

DO 1C 1 = 1 , 1C 
DO 10 J = 1 , 2 0 

T ( I , J ) = 0. 

10 R ( I , J ) = C. 

C 


FOC 

- 

X2Y 



F 1C 

- 

X2Y?*< 1 * 

+ 

SE 1) /fc. 

FG1 

= 

X ?Y?/6. . 



F2(> 


X2Y?* ( 1 . 

-H 

5 E 2 S fc 2 ) /I 2 . 

FI J 

- 

X 2 V 3=M 1 . 


2 . * S E 1 5 /2a • 

FO? 


X ?Y3/1 ? • 



F30 

- 

X 2Y?* ( 1 . 


SE 2 •» St 2 S.E3)/2C • 

F 2 1 

- 

X 2Y ?.* ( 2 . 

+ 

?.* SE1 + 3 ,*SE2 )/60. 

FI? 


X2Y?*(1 . 

♦ 

. 3 »*=S E 1 )/60. 

FC? 


X?Y3/20 . 



FAC 


x?v:<*n . 

+ 

S t 1 SE2 ♦ SI 2 ♦ SFa)/30. 

FrI 


X 2 Y3* ( 1 . 


?.*sn ♦ 3.*sr2 ♦ A.*sr?)/i 

F2? 


X?Y3*U . 


3.*SE1 fc.*SE2)/lF0. 
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FI? = X2Y3*Cl. -A-.*S£1 )/120* 


FOA = X2Y3/30. 

Tf ltl» 

= 

FCC 

TCI,?) 

= 

Fie 

m.3i 

- 

FC1 

TCI, a) 


FI 1 

TC l*f 1 


F02 

T C 2,2 ) 

“ 

F 20 

T C 2 * 3 ) 

= 

F 1 1 

TC?,a) 

r 

F ?! 

T C 2 , 5 ) 

= 

F 1 2 

TC? ,3 > 

JE 

FO? 

T(?,a) 

= 

FI? 

T C 3 * F ) 

- 

FC? 

7Ca,a) 

= 

F2? 

TCa, 5) 

= 

F 1 3 

TCF,f ) 

=• 

FC a 

TCfr ,fc) 

- 

F( 0 

T C fc ,7 1 

TT 

FIG 

TC6,P) 

- 

FC1 

TCfc ,*») 

- 

F?C 

T C 6 , 1 0 ) 


= F 1 1 

T C 7 ,7) 

= 

F2C 

T C 7 , 6 ) 

- 

FI! 

T C 7 , ° ) 

=: 

F3C 

TC7,K ) 


= F?1 

TC £ ,F ) 

- 

Ft? 

TCP, 0 ) 


F ? 1 

TCF.IC) 


= F 1 2 

TC°,o) 


FAC 

T C ° , 7 0 ) 


= F?7 

T CIO, 1C ) 

= F22 

DC 30 1 


7,IC 

DO 30 J 


I .10 

TCJrl) 

=■ 

T 1 1, J 


= 1. 

= -1. 

= 1 . - 
= sri - i 
= - v ? 

= -St! 

- 7 . 

= >3 
= Y2 

- — Y2 

- Y?*C 1 . 

= Y?~St 1 

=■ -Y? 

= 1. 

= —7 . 

= >2 
-- 1 * 

= -Y? 

= S.L 1 - 


7 »7I 
2 , 1 ) 

2 ,aj 
3,1) 
3,?) 

?,7> 

3,°) 

4 . 3 ) 
A ,6) 

3.3) 
c ,6 ) 

3 ,°) 

f>»? ) 
•>,2) 

?,3) 

7.3) 

7,h) 

fc ,2 ) 


SF1 ) 


1 



H2A2 


3/ 3 


M P * 5 ) = -It l 
RtP.fel = X? 

R(P,P) = 1. 

F (P *° ) = -X? 

R ( ^ »3 ) = -X2 
F(°,6) = x: 

P ! 1 0*2 ) = X2*(SE1 - 1.1 
P(lu.6) = -X3 
-R ( 10*9 ) = X.2 
C 

CAUL PTAEA CT,R,10,9,ja ,KR1 . 
DC -0 I-l,° 

nr j=i,° 

40 ZII.J1 = TCI » J ) 

C 

- RETURN 
END 
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SUBROUTINE M2B1 <X2 ,Y3, TH,RG,Z,KZ I 
DIMENSION Z(KZ,1) 

C 

C SUBROUTINE Tf CALCULATE FINITE ELEMENT... 

C LUMPED MASS MATRIX* 

C FPR A eFNDING TRIANGLE PLATE ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C MASS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE PLATE TO LIE IN AN X-Y PLANE 
C WITH JOINT I AT THE X-Y CPIGIN, JOINT 2 LIES ALONG THE POSITIVE 
C X AXIS, AND JOINT 3 IS IN THE POSITIVE Y DIRECTION. 

C LOCAL COORDINATE ORDER IS \ 

C CDZ,TX,TY ) JOINT 1, THEN JOINT 2, 3. 

C WHERF DZ IS TRANSLATION AND TX,TY ARE ROTATIONS.^ 

C DEVELOPED FY WA B ONFIELD • FEBRUARY 1973. ~ > 

C LAST REVISION BY RL WOHLEN. SEPTEMBER 1973. 

C 

C SUBROUTINE ARGUMENTS 


c 

X2 

= INPUT 

LOCAL X COORDINATE LOCATION OF 

JOINT 2. 

C 

Y3 

= INPUT 

LOCAL Y COORDINATE LOCATION OF 

JOINT 3. 

C 

TF 

= INPUT 

PLATE THICKNESS. 

• 

c 

RO 

= INPUT 

MASS DENSITY. 


C 

Z 

= OUTPUT 

MASS MATRIX. S1ZE(9,9|. 


c 

KZ 

= INPUT 

ROW DIMENSION OF Z IN CALLING 

PROGRAM. 


C 

AREA = C.B*X2*Y3 
CM = (RC*TF*AREA)/3.0 
DO 10 1=1 , c 
DC 1C J=l,o 
10 2(1, J) = 0.0 
DC 20 1 = 1 ,* 

20 ZII,1) = CM 
RETURN - 
END 
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SUBP.CUT INF M2F2 ( X? ,X 3, Y3 ,TH,FHC, Z,T,P ,KZ ,KT,KR ) 

DIMENSION 2 (KZ *11*1 (KT * 1 ) ,R (KR ,11 " 

DIMFNSICN IVEC ( 1C ) »COEF 19 ) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C CONSISTENT NASS MATRIX, 

C FOR A FENDING TRIANGLE PLATE FLEMENT WITH UNRESTRAINED BOUNDARIES. 

C CUBIC DISPLACEMENT (LINEAR CURVATURE ) FIELD IS USED. 

C THIS IS NOT ThE St* CALLED STRICKLAND ELEMENT. 

C MASS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C ThE LOCAL COORDINATE SYSTEM ASSUMES THE PLATE TO LIE IN AN X-Y PLANE 
C WITH JOINT I AT THE X-Y ORIGIN* JOINT 2 LIES ALONG THE POSITIVE 
C X AXIS, AND JOINT 3 IS IN THE POSITIVE Y DIRECTION. . 

C LOCAL COO c DINaTE CFDFR IS 
C (DZ,TX,TV) JOINT 1, THEN JOINT 2, 3. 

C WHERF OZ IS TRANSLATION aND TX,TY ARE ROTATIONS. 

C CALLS FORMA Sl*F°OL‘TlNES BTABA. 

C DEVELOPED EY CaF L FODLEY MARCH 1«75. 

C LAST REVISION EY WA BENFIELD. SEPTEMBER 1°?3. 

C c 

C SUBROUTINE ARGUMENTS 


c 

X2 

=1 

I^PUT 

LOCAL X COORDINATE 

LOCATION OF 

JOINT 

2. 


c 

X3 


INPUT 

LOCAL X COORDINATE 

LOCATION OF 

JOINT 

3. 


c 

Y3 


lttPtrr 

LOCAL Y CCCPOINATE 

LOCATION OF 

JOINT 

3. 

: 

c 

TH 


INPUT 

PLATE THICKNESS. 





c 

RHC 


INPUT 

MASS DENSITY. 





c 

Z 

= 

CVTPV1 

MASS MATRIX. SIZE(9,9). 




c 

T 

= 

INPUT 

MATRIX WORK SPACE. 

SIZF(10,10) 

• 



c 

P 

- 

INPUT 

MATRIX WORK SPACE. 

SIZF ( 10,101 

• 



c 

KZ 


INPUT 

Ft* DIMENSION OF Z 

IN CALLING 

PROGR AM. 

MIN=9. 

C 

KT 

= 

INPUT 

ROW DIMENSION (F T 

IN CALLING 

PPCGR AM. 

MIN= 10. 

c 

KR 


INPUT 

ROW DIMENSION OF R 

IN CALLING 

PROGRAM. 

MIN=!0. 


C 

ssi a x?/x? 

SS 2 = SF 1 *SF 1 
SE? = SE 2 *SE I 
SE^ = 5 F 3 *SE I 
SEE = SENSEI 
SFfc = SEE*SE 1 
SEC 1 = Cl. + StlI/ 3 . 

srci = s c ci**? 

SECS = SEC 1**2 

c 

ro ic i = i ,k 

DO 1 C J= 1 , 1 C 
10 T (I ,J) = 0. 

C 

T (1 , 1 ) = I . 

T ( 2 , ? ) = 1 . 

T ( ? , 2 ) = i. 

T(A, 1 ) = 1 . 

1 U, 2 ) = 1 • 

T ( a t a ) = 1 . 

T ( A ,7 ) = 1 . 

T ( E » 3 ) = 1 . 
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T( 5 » 5 1 


l. 

7 ( 5 »P ) 

- 

i. 

7(6,2) 

= 

i. 

7(6,4) 

- 

2 • 

7 ( 6 ,7 ) 


3. 

7(7,1) 

=■ 

1. 

7(7,2) 

= 

SB 1 

7(7,3) 

- 

1 . 

7(7, M 


SE2 

7(7,5) 


SE1 

7(7,6) 

- 

1. 

7(7,7) 

=: 

SF? 

7(7,6) 


SF? 

7(7,9) 


SF1 

7(7, 1CI 

1 5 

= 1. 

7(8,31 


1 . 

7(6,5) 

= 

SE1 

7(8,6) 

~ 

2. 

7(8,6) 


SE? 

7(P,oi 

— 

2.*SE1 

T (8 , 10 ) 

= 3. 

7(9,2) 

— 

1. 

7(9,4) 

= 

?.*SE1 

7(0,5) 

=r 

1 . 

7(0,7) 

- 

?.*SE? 

7(9,8) 


2.*SE1 

T( 9, 9 ) 

= 

1. 


7(10,1) = 1. 

T ( 10,2 ) = SEC1 
7(10,3) = 1./3. 

7(70,4) = 5FC2 
7(10,5) = SEC1/2*. 
T|10,6) = 1./9. 

7(10,7) = SEC? 

T ( 1 0 , f ) = SEC2/3. 

7(10, o) = SEC 1/9 • 
7(10,10) = 1./27. 

C 

DC £ I = 1»H* 

DC 7 J=1 ,10 
7 R f I ,J) = 0. 

5 R(I,I) = 1. 

C 

DC 1 GO L-l ,10 
JBIG =• 1 

A 1 =• ABS(T(L»1)) 

DO 15 J = 2 ♦ 1 0 
A? = A£S(7(L,J)> 

IF ( A2 .LT. Al) GO TO 15 
A1 = A? 

JBIG = J 
15 CONTINUE 

IVFC(l) = JBIG 
ALJEIG = Til, JBIG) 

DC 17 J=1,1C 
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TCL,J) = T(L,J)/ALJFIG 
17 R(L,J) = F (L,J)/ALJB1G 

nr 25 1=1,10 
AIJFIG = TdfJUGI 
IF (I • FQ . L) GC TC 25 
DO 30 J=1 *10 

T ( I , J 1 = HI,J> - AIJB1G*T(L, J) 

30 R ( I » J 1 = RCI,JJ - AIJBIG*R(L, J) 

25 CONTI NUF 
100 CONTINUE 

DO 1*1, 1C 
IP = IVFCII) 

DC *0 J = 1 »!0 
40 T ( 1R » J 1 = P(I,J) 

DC 50 1 = 1 ,10 
rc 5o j=i,it 
50 R(I,J) = T II , J ) 

DO 20 1=1,1 ( 

R ( I » 2 1 = Y?*P(I,21 

PC 1,31 = — X2*F (1,3) 

R ( I *5 ) = Y3*P.(I,5) 

PCI, 6) = — >'2*R C I ,6) 

R C I , £ ) = Y3*F(I,8) 

20 R ( I ,9) = -X2*M1»V|- 

COEF(l) = 1./3. 

C0EFC2) = Y3/1P. 

COFFC?) = -CX2+X31/16. 

C0EFC4) = 1./3. 

C0FFC5) = Y3/1E. 

COFFC 6 1 = C2.*X2 - X3 1/16. 

C0EFC7) = 1./3. 

CCFFC8) = -Y3/«. 

C0EFC9) = C2.*X3 - X2 1/1B. 

DO 80 1 = 1 ,10 
DO 80 J=l,9 

80 K C 1 , J) = R ( 1 , J ) + RCI ,10)=»COEFCJl 

DO 55 1=1,10 
DO 55 J=1,!C 
55 TCI, J) = 0. 


X2Y3 


= X?*Y?*7K*RH0 

FGo 


X2V3/2. 



FID 


X2Y.VM1 . 

4 

SE 1) /6. 

PCI 


X ?Y 5/6 . 



F?0 

=: 

X?Y r =*(l . 

4 

SE1 + SE21/12. 

Ell 

2: 

X2Y;>*(1 . 

4 

2 .*SC1 )/24. 

F02 

- 

XPY'Vl? . 



F 30 

- 

x?v?*ci . 

4 

i t 1 4 S E 2 + f E 3 1/20 

F21 


X2Y 3* ( 1 • 

4 

2 .*S FI + 5.*SE2)/60 

FI? 


X2Y3* tl . 

4 

3.*5E1 1/60. 

F03 


X 2 Y 3/20 * 
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F4C = X2Y3M1. + SE1 ♦ SE2 ♦ SE3 ♦ SE4)/30. 

F31 = X2Y3*(1. ♦ 2.*SF1 3.*SE2 ♦ 4.*SE3)/I20. 

F22 = X?Y3*(1. + 3.*SE1 6.*SE2)/180. 

F 13 = X2Y3*<1. + 4.*SE1)/12G. , 

F04 = X2Y3/30. 

F. r 0 = X?Y3*fl. + SE1 ♦ SE 2 + SF3 ♦ SE4 ♦ SE51/42. 

F41 a X2Y3* ( 1 . ♦ 2 »*S El ♦ 3.*SE2 + 4.*SF3 ♦ 5.*SE*)/210. 

F32 = X?Y3*(1. + 3.*SE1 + 6.*SE2 + 10.*SE3l/420. 

F23 = X2Y3*(1. ♦ 4.*SE1 ♦ 10.*SE2 )/420. 

F 14 = X?Y2*Cl. + 5.*SE1)/21G. 

F 05 = X2Y3/*>2. 

F60 = X2Y3*(1. SE 1 ♦ SE2 ♦ SE3 ♦ SE4 + SE5 + SE6)/56. 
F51 = X2Y3*(1 .*2.*SE) +3.*SE2+4.*SE3«-5.*SE4+6.*SE5 )/336. 
F42 * X2Y3 S M1 •+3.*SE1 ♦6«*SE2-*‘1C .*SE3-*T5 .*SF4) /840 • 

F33 = X2Y 3* 1 1 . +4. *SE1 +1 0* * : SE2 *20 »*SE3 1/1120 • 

F24 = X2Y?*(l.+E.*SEl+lb.*SE2)/840. 

F15 = X2Y3*I1 .+6.*SE1 1/336. 

F06 = X2Y2/56. 

C o 


Til ,11 

- 

F00 

T<1, 2) 

= 

FI C 

T(1 ,3) 

= 

F01 

T(1 ,4) 

- 

F2C 

T( 1 * 5 ) 


F 1 1 

Ttl v 6) 


F02 

TCI, 7) 

=■ 

F3C. 

TCI, 8) 

- 

F?1 

T(l,9) 

- 

FI 2 

T(1 ,10) 


F03 

T( 2 ,2 ) 

rr 

F2C 

T ( 2 ,3 ) 

= 

FI 1 

T( 2 ,4) 

— 

F30 

T ( 2 , f> ) 


F21 

T C 2 *fc ) 

— 

F12 

T ( 2 » 7 ) 

— 

F40 

T (2,8) 

- 

F31 

T ( 2 , 9 ) 

= 

F22 

T C 2 , 1 0 ) 

= 

F13 

T( 3 , 3 ) 


FC'2 

T ( 3 ,4 ) 

- 

F21 

T C 3 , 5 ) 

- 

F 1 2 

T ( 3 ,6 ) 

- 

F03 

T ( 3 ,7 ) 

=r 

F?1 

T ( 3 , 8 ) 

=■ 

F22 

T(3,o) 


FI 3 

T(3 ,10) 

= 

: F 04 

T ( 4 , A- ) 

=. 

F40 

T (4 , 5 ) 


F3 ) 

T(4,6) 


F22 

T (4,7 ) 

- 

F SO 

T(4,E) 


F A 1 

T ( 4 , 9 ) 

- 

r?? 

T ( 4 , 1 (! ) 

= 

F ?3 

T(F,5 ) 


F22 

T(i>,o) 

- 

F13 
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C 


C 


T ( 5 ,7 ) 

= 

F^l 

T (5,8) 


F32 

1(5,9) 


F23 

T ( 5 » 1 0 ) 


= F 14 

T ( 6 ,6 ) 


F04 

T (6,7) 


F?2 

T ( 6 , 6 ) 

= 

F?3 

T (6,9) 

n 

F14 

T (6,10) 


= F05 

T ( 7 ,7 ) 

- 

FfcC 

T(7,f ) 

=r 

F51 

T ( 7 ,9 ) 


F^2 

T ( 7 , 10 ) 


= F33 


T ( R ,8 ) = F^? 

T(6,9) = F33 
T ( 8 , 1 0 ) = F24 
T(<?,°) - F?4 

T ( ° , 1 0 ) = F 1 5 
T( 10,10) .= FCfc 

DC 60 1=1,10 N 
DC 60 J=I,10 
60 T ( J , I ) = T ( I » J ) 

CALL HTAPA CT,R,10,9,KT,KR) 
DO E5 1 = 1 ,9 , .-= ' 

DO 85 J= 1 ,9 
65 1(1, J) = 1(1, J) 

RETURN 

END 



M3C1 


Silt PONTINE M3C 1 (X3,Y3 ,TH,RO ,2 ,KZ ) 

DIMENSION Z ( KZ » 1 ) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C LUMPED MATS MATRIX 

C FOR A. RECTANGULAR 5 HF AR PANEL ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C MASS MATRIX IS IN LOCAL COORDINATE SYSTEM. 

C THE LOCAL COORDINATE SYSTEM ASSUMES THE PANEL TO LIE IN AN X-Y PLANE 
C WITH JOINT 1 AT THE X-Y ORIGIN, JOINT 2 LIFS ALONG THE POSITIVE 
C X AXIS, JOINT ? IS IN THF POSITIVE X,Y DIRECTION, AND JOINT 4 LIES 
C ALONG THF POSITIVE Y AXIS. 

C LOCAL COORDINATE ORDER IS 
C DXl,DX2,DX3,DX* f DY 1 , DY 2 , LY3 * DY4 

C WHERE DX,DY ARE TRANSLATIONS. 

C DEVELOPED BY RL WOHLEN. APRIL 1974. 

C 

C SUBROUTINE ARGUMENTS 

C X3 = INPD. LOCAL X COORDINATE LOCATION OF JOINT 3. 

C Y? = INPUT LOCAL Y COORDINATE LOCATION OF JCINT 3. 

C TH = INPUT PANEL THICKNESS. 

C RO = INPUT MASS. DENSITY. 

C Z = OUTPUT MASS MATRIX. SIZE(&,8). 

C KZ = INPUT ROW DIMENSION .OF Z IN CALLING PROGRAM. MIN-8,. 

C 

CM - RO*TH*X2*Y3/4.0 
DO 10 J = ! ,E 
DO 10 1=1 ,t 
10 Z(I,J) = O.C 
DO 20 1=1, E 
20 2(1,1) = CM 
C 

RETURN 

EN'' 
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SUBROUTINE MAS1A t CJ ,E J , Al ,A2,RU , ‘JAMEM ,Z *KC J ,KEJ ,KZ ) 
DIMENSION CJ(KCJ.l), £J(KEJ,1>, 2 : . ' * 1 > 

DIMENSION E 1 i 3 *3 1 * E2(3,3I, W(2,2) 

/ 

SUBROUTINE TO CALCULATE FINITE ELFMENT... 

MASS MATRIX 

FOR AN AXIAL ROD ELEMENT WITH UNRESTRAINED BOUNDARIES. 
ROD MAY BE LINEARLY TAPERED OR UNIFORM. 

MASS MATRIX IS IN GLOBAL COORDINATE DIRECTIONS. 

GLOBAL COORDINATE ORDER IS 

(U,V,W) JOINT 1, THEN JOINT 2. 

WHERE U,V,W ARE TRANSLATIONS. 

EULER ANGLE CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 
CALLS FORMA SUBROUTINES EULER, Ml Al ,M1 A2 * ZZEOMB . 

DEVELOPED BY RL WCHLEN. SEPTEMBER 1972. 

LAST REVISION BY WA BENFIELD. MARCH 1976. 

SUBROUTINE ARGUMENTS 


CJ 


INPUT 

MATRIX OF GLOBAL X,Y,Z COORDINATES AT ROD JOINTS. 
POWS 1,2,3 CORRESPOND TO X,Y,Z COORDINATES. 

COLS 1,2 CORRESPOND TO JOINTS 1,2. SIZE (3,2). 

EJ 


INPUT 

MATRIX OF EULER ANGLES (DEGREES) AT ROD JOINTS. 
ROWS 1,2,3 CORRESPOND fO GLOBAL X,Y,Z PERMUTATION. 
COLS 1,2 CORRESPOND TO JOINTS 1,2. SIZE(3,2). 

Al 

= 

INPUT 

CROSS-SECTION AR FA AT ROD END 1. 

A2 


INPUT 

CPOSS— SE CTION AREA AT POD END 2. 

RO 

= 

INPUT 

MASS DENSITY. 

NAMEM 


INPUT 

TYPE OF MASS MATRIX WANTED. 
= Ml, LUMPED. 

= M2, CONSISTENT. 

Z 


OUTPUT 

MASS MATRIX. SIZE(6,6). 

KCJ 

= 

INPUT 

ROW DIMENSION OF CJ IN CALLING PROGRAM. MIN=3. 

KEJ 

= 

INPUT 

ROW DIMENSION OF EJ IN CALLING PROGRAM. MIN=3. 

KZ 

= 

INPUT 

ROW DIME NS ION OF 1 IN CALLING PROGRAM. MIN=6. 


NERROR EXPLANATION 

1 = DIMENSION SIZE EXCEFDED (KZ) . 

2 = NAMEM IMPROPERLY DEFINED. 


IF (KZ .LT. 6) GO TO 999 
DO 6 J=l,6 
DO 5 1=1,6 
5 Z(I,J) = 0.0 

RL = SQRTUCJ(1,2)-LJ (1,1 >>**2 ♦ ( C J ( 2, ? >-C J ( 2 , 1 ))**2 
* + (CJ<3,2>-CJ(3,\))**2) 

IF (NAMEM ,E0. 6HM1 ) GC TC' 110 

IF (NAMEM .FQ. 6HM2 ) 00 TC 120 

GO TO 999 


NERR0R=1 


NERR0R=2 


LUMPED. 

110 CALL Ml A 1 ( Al,A2,RL,RC,W,?) 

DC 112 1*1,3 
112 Zf 1,1 ) = W(l,l 1 



n o 
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DO 114 1=4,6 

114 Z(I,n = W(2, 2 ) - 
RETURN 

CONSISTENT. 

120 CALL Ml A2 < A1,A2,RL ,R0,W,2) 

DO 122 T = 1 » 3 
122 Z(I,I) = Wtl.l) 

CALL EULER ( E J C 1 , 1 ) , El ,3 ) 

CALL EULER < E J < 1 ,2 > » E2 ,3 ) 

CALL ATXBB (E1,E2, 3,3,3, 3,3) 
DO 124 1=1,3 
DO 124 J=4,6 

ZtI,J) = W(1,2)*E2(I» J-3) 

124 ZIJ,I> = ZII,J) 

DO 126 1=4,6 
126 Z(I,I) = W(2,2 ) 

RETURN 

999 CALL TZBL’MB .-( 6HMAS1A ,NERROR) 
END 



MAS IB 


1 / 2 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

r 

V 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE HAS1B (CJ,EJ M,A2 ,R I 1 ,P 12 ,RO,NAME M, Z, W,KC J ,KB J, K2 »KW ) 
DIMENSION C J( KCJ » I ) , EJ(KEJ,1), 2 CKZ, 1) , W(KW* 1) 

SUBROUTINE T D CALCULATE FINITE ELEMENT... 

HASS MATRIX 

FUR A COMBINED AXIAL— TORS ION-BENDING BAR ELEMENT WITH UNRESTRAINED 
BOUNDARIES. 

BAR MAY BE LINEARLY TAPERED UR UNIFORM. 

MASS MATRIX IS IN GLOBAL COORDINATE DIRECTIONS. 

GLOBAL COORDINATE ORDER 11 

(U,V,W,P,C,R) JOINT 1, THEN JOINT 2 
WHERE U,V»W ARE TRANSLATIONS AND P,Q,R ARE ROTATIONS. 

EULER ANGLE CONVENTION IS GLOBAL X,Y,Z PEP MUTATION • 

CALLS FORMA SUBROUTINES BT AB A, DCOS1B , M1A1, Ml A2 ,M1B I* M1B2 ,M1C 1 ,MIC2 , 

ZZPOMB . 

DEVELOPED BY RL WOHLEN. FEBRUARY 1973. 

LAST REVISION BY WA BENFIELD. MARCH 1976. 

SUB ROUT INF ARGUMENTS 

INPUT MATRIX OF GLOBAL X,Y,Z COORDINATES AT BAR JOINTS. 

ROWS 1,2,3 CORRESPOND TO V ,Y,Z COORDINATES. 

COLS 1,2 CORRESPOND TO oOINTS 1,2. COL 3 CORRESPONDS 
70 REFERENCE POINT TO DEFINE LOCAL XY PLANE. SIZE{3,3). 
INPUT MAT^ty CF EULER ANGLES (DEGREES) AT BAR o. iNTS. 

ROWS 1,2,3 CORRESPOND TO GLOBAL X,Y, Z PE RE UTATION. 

COlS 1,2 CORRESPOND TO JOINTS 1,2. SIZE(3,2). 

INPUT CROSS-SECTION AREA AT BAR END 1. 

INPUT CROSS-SECTION AREA AT BAR END 2. 

INPUT CFOS S-SE CT ION POLAR AREA MOMENT OF INERTIA AT END 1. 

INPUT CROSS-SECTION POLAR AREA MOMENT CF INERTIA AT END 2. 

INPUT MASS DENSITY. 

INPUT TYPE CF MASS MATRIX WANTED. 

= Ml, LUMPED. 

= M2, CONSISTENT. 

OUTPUT MASS MATRIX. SIZE (12,12). 

INPUT WORK SPACE MATRIX. SIZE(12,I2). 

INPUT PCW DIMENSION OF CJ IN CALLING PROGRAM. MIN=3. 

INPUT ROW DIMENSION C r EJ IN CALLING PROGRAM. MIN=3. 

INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=12. 

INPUT ROW DIMENSION OF W IN CALLING PROGRAM. MIN=12. 

N'FRROP EXPLANATION 
1 = DIMFNSICN SIZE EXCEEDED KZ,KW). 

?. = NAMEM IMPROPERLY DEFINED. 


CJ 


EJ 


A1 

A? 

PI1 
P 12 
RO 

NAMEM - 


Z 

W 

KCJ 

KEJ 

KZ 

KW 


NERR0R=1 

IE (KZ . LT 4 12 .OP. KW .LT. 12) CO TO 999 
DC f- J = 1 , 1 2 
DO 5 1=7,12 
5 Z(I,J) = c.O 

RL = SORT (CJ (I ,2 )-( J (1 ,1 ))**2 ♦ ( C J ( 2, 2 )-C J ( 2 , 1 ) )**2 
* ♦ (CJ(3,2)-CJ(3,1 ) )**2 ) 

IF (NAMEM ,Fv. . OHM 1 > GO TO 11C 

IF (NAMEM . PC • 6HM2 ) EC TC 120 


NERROR "2 



o n 
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GO TO ®°5 


C 

c 

c 


e 


AXIAL=«1A1 


11C 


Call mIaI 
•- ALL M1C1 
CALL *151 
CALL MJR1 

GO TO 3C0 



bend:\!&=hiei 


AXIAL =W1A2 (LINEAR DISP), T0RSI0N=M1C2 I LINEAR DISP), 
BE\*DING=M1F? CCUFIC DISP ). 

120 CALL m?A 2 f A 1 , A2 ,RL ,PC , L,KZ ) 

CALL M1C2 (P.U,PI2,RL,RC,Z(3»3),KZ) 

CALL MP2 f A1,A2,RL,P0,ZI5,5),KZ1 
DO 125 J=7»6 
DO 125 I=5 t 6 
ZCI.J) =-ZII,J) 

125 Z(J,I) =— Z (J, I ) 

CALL M1F2 ( A1,A2,PL ,P0,Z(9,9),KZ) 

300 CALL DCCS1F (CJ,EJ,K,KC J,K.F J,KW) 

CALL ETAPA (Z,W, 12,12, KZ,KW) 

RETURN 


! LUMPED )® 


C 


999 CALL ZZFC.Mf (6HP.AS1E ,NERPPP) 
END 



oonooooooooononooooooriooonn 


MAS2 


1 / 2 


SUBROUTINE MAS 2 ICJ,EJ,TMAS ,RC,NAMEM,Z,Wl ,W2,KC J,KE J,KZ ,KW1,KW2 I 

DIMENSION CJ(KCJ,1), EJ(KEJ,1), Z(KZ,1>, W1(KW1,1), W2(KW2,1) 

DIMENSION IVFCC1FI 

DATA 1VFC/I,2»6,7,8,12,13*l^,ie, 3*4,5,9,10,11,15*16*17/ 

SUBROUTINE TC CALCULATE FINITE ELEMENT... 

MASS MATRIX 

FOR A COMBINED MEMBR AN E- FENDING TRIANGLE PLATE ELEMENT WITH 
UNRESTRAINED BOUNDARIES. 

MASS MATRIX IS IN GLCFAL COORDINATE DIRECTIONS. 

GLOBAL COORDINATE ORDER IS 

(U,V,W,P,t,R) JOINT 1, THEN JOINT 2* 3. 

WHERE U,V,W ARE TRANSITIONS AND F,0*R ARE ROTATIONS. 

EULcR ANGLE CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

CALLS FORMA SUBROUTINES pTAEA,DC 0S2,M2A1 ,M2A2,M2B1 ,M2e2,ZZE0MB . 
DEVELOPED EY WA FENEIELD » RL WGhLEN. FEBRUARY 1973. 

LAST REVISION EY Vs BENFIFLD. MARCH 1976. 


C 

C 

C 

c 

c 

c 

c 

c 

c 

c 


SUE ROUTINE 


CJ 

• — 

INPUT 

EJ 


INPUT 

TMAS 

_ 

INPUT 

RO 


IN^UT 

NAMEM. 

- •; 

INPUT 

Z 

_ 

OUTPUT 

W1 

- 

INPUT 

w? 

- 

INPUT 

KCJ 

=: 

INPUT 

KEJ 

- 

INPUT 

KZ 

= 

INPUT 

KW 1 

- 

INPUT 

KW2 

= 

INPUT 


ARGUMENTS 

MATRIX OF GLOBAL X,Y,Z COORDINATES AT TRIANGLE JOINTS. 
ROWS I ,2 ,3 CORRESPOND TO X,Y,Z COORDINATES. 

CCLS 1,2,3 CORRESPOND TO JOINTS 1,2,3. SIZEC3.3I. 
MATRIX OF EULER ANGLES (DEGPEFSI AT TRIANoLE JOINTS. 
ROWS 1,2,3 CORRESPOND TC GLOBAL X,Y,Z PERMUTATION. 

CCLS 1,2,2 CORRESPOND TO JOINTS 1,2,3. SIZE(3,3J. 
EFFECTIVE MACS THICKNESS. 

MA^S DENSITY. 

TYPE OF MaSS MATRIX WANTED. 

= Ml, LUMPED. 

= M2, CONSISTENT. 

MASS MATRIX. SIZEC1B,181. 

WORK S FA 05 MATRIX. SIZF{1B,16). 

WORK SPACE MATRIX. $IZE(J0,10). 

ROW DIMENSION OF CJ IN CALLING PROGRAM. M1N=3. 

ROW DIMENSION OF EJ IN CALLING PROGRAM. MIN=3. 

COW DIMENSION OF Z IN CALLING PROGRAM. MIN=18. 

ROW DIMENSION OF W1 IN CALLING PROGRAM. MIN=1B. 

COW DIMENSION OF W2 IN CALLING PROGRAM. MIN=10. 


NE C FTP EXPLANATION 

1 = D IH C NS IC'N SIZE EXCEEDED (K 2 , KW* .KV2 ) . 

2 = NAME* IMPROPERLY DEFINED. 


NEPR0R=1 

IF <KZ .L7. IE .OR. KW! .LI. lb .OR. KW2 .LT. 10) GO TO 999 
DC 5 J=l,lb 
DC E Jrlt'P 
5 Z C I * J ) - < .C 

S L 1 2 = 5C C H(CJ<1,2)-CJ(1,1))**2 ♦ <C J( 2 , 2 )-C J ( 2 , 1 )) **2 

* ♦ (CJ(?,?)-CJ(3,l))*-*2) 

SL22 = SC PTKCjn ,? >— C J 1 1 ,2 ))**2 + <C J( 2,3 ) -C J ( 2, 2 ) ) **2 

* ♦ (CJ(3,3)-CJt3, ?))**?) 

S L 13 - f C 1 T ( ( C J ( 1 ,2 )-CJ 1 1 , 1 ) ) **2 ♦ If Jf ?» 3)- C J C 2, 1 ) )**2 

* ♦ (CJ(3*3)— CJ(3»1) 1**2 ) 

X3 = <SL13’?*2 + SH2**2-SL23**2)/<2.0*SL12I 



no on 


HAS? 
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Y3 = SC>STISLT3**2-X3**2) 

IF (NAMBM .EC. fcHHl J GT TC 110 

IF (NAMEM .EC. 6HM2 ) GC TO 120 

GO TO 999 

MEMBRANE = M 2A1 I LUMPED I , FENDING = M2F1 (LUMPED). 
110 CALL M2 Al C ? L : 2 ,Y3. TMAS ,F 0, Ml ,KK1 ) 

CALL M2F1 ( SL12,Y3,TKAS,F'C,W1(10,10),I‘W1) 

DO 115 IW = 1 * 1 E 
IZ = IVFCIIK) 

115 ZfIZtIZI = WltlWtlW) 

RETURN 


NERR0R=2 


MEMBRANE = M?A2 (CONSISTENT), BENDING = M2E2 (CONSISTENT). 

120 CALL M2A2 ( SL12 ,> 3, Y3 ,TMAS, RO,Z ,W1 ,W2 ,K Z.KWl ,KW2 ) 

CALL M2F? ( SL12 ,X3, Y? ,TMAS ,PD,Z ( 10, 10 ) ,Wl ,7f 2,KZ *KW1 ,KM? ) 

CALL JC0S2 (CJ,EJ,W1 , KCJ,KFJ ,KH1 ) 

CALL FTABA (Z,Wl ,18, If ,KZ,KW1) 

RETURN 

999 CALL ZZEPMB (5H.HAS2 ,NER?OR ) 

END 
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SUP ROUT IN r MAS3 (CJ,E J,TMAS,R0,NAMEM,Z,W1, 2,KCJ,KEJ,KZ,KW1,KW2) 

DIMENS I UN CJ(KCJ,I»,E J ( KE J, 1 ) , Z (KZ , 1 ) ,W1 (KWll l ) ,W2(KW2,1 ) 

DIMENSION CW(3,2), EW(3,3), 1V1(18), IV2(18)$, IV3(ie>, IV4(18), 

* W3( 1C *1C) 

DATA IV1/ I* 2, 3, 4, S, t>, 7, 6* ? ,1 C t ! 1 *1 2, 13,1*, 15. 16*17*18/ , 

* IV2/ 1* 2, 3, <*, 5, 6,13,14,15,16,17,18,19,20,21,22,23,24/, 

* 1V3/ It 2, 3, 4, 5, 6, 7, B, 9,10, 11 ,12, 19,20,21 ,22,23,24/, 

* 1V^/ 7, F , 9,10, 11 ,12, 13,14,15,16,17,18,19,20,21,22,23,24/ 

C 

C SUBROUTINE T P CALCULATE FINITE ELEMENT... 

C MASS MATRIX 

C FOR A COMBINED MEMBRANE-PENDING QUADRILATERAL P L ATE ELEMENT WITH 
C UNRESTRAINED E OUND ARIES. 

C MASS MATRIX IS IN GLOBAL COORDINATE DIRECT IONS . 

C GLOBAL COORDINATE ORDER IS 

C (U,V,W,P,C,R) JOINT 1, THEN JOINT 2, 3, A. 

C WHERE U,V,W ARE TRANSLATIONS AND P,C,R ARE DOTATIONS. 

C EULER ANGl r CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

C CALLS FORMA SUBROUTINES MAS2 ,EFVADD,ZZPOMP . 

C DEVELOPED B Y WA F.ENF IELC , RL WCHLEN. FEBRUARY 1973. 

C LAST REVISION BY WA BENE IE LD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS 

C CJ = INPUT MATRIX OB GLOBAL X,Y,Z COORDINATES AT QUAD JOINTS. 

C ROWS 1,2,3 CORRESPOND TO X,Y ,7 COORDINATES. 

C CCLS 1,2 ,3,4 CORRESPOND TO JOINTS 1,2, 3, 4. SIZEI3.41. 

C EJ = INPUT MATRIX OF EULER ANGLES (DEGREE SI AT QUAD JOINTS. 

C ROWS 1,2,3 CORRESPOND TO GLOBAL X,Y,Z PEPMUTATION. 

C COLS 1,2 ,3,4 CORRESPOND TO JOINTS 1, 2,3,4. SIZE(3,4). 

C TMAS = INPUT EFFECTIVE MASS THICKNESS. 

C RC = INPUT MASS DENSITY. 

C NAMEM = INPUT TYPE CF MASS MATRIX WANTED. 

C = Ml, LUMPED. 4 TR I ANGLES, OVERLAP AVERAGE. 

C = M2, CONSISTENT. 4 TRIANGLES, CVFRLAP AVERAGE. 


C 

■» 

d. 

=: 

UTPUT 

MASS 

MATRIX. S 

IZE (24,24) . 



C 

W! 

= 

I VP UT 

WORK 

SPACE MATRIX. 

SIZF(18,18) 



C 

W2 

2 

INPUT 

WORK 

SPACE MATRIX. 

SIZE (18,18) 

• 


C 

KCJ 


INPUT 

ROW 

DIME NS ION 

OF CJ 

IN CALLING 

PROGRAM 

. MIN=3. 

c 

KE J 

= 

INPUT 

FPW 

DIMENSION 

OF FJ 

IN CALLING 

PRCO-R AM 

. M1N=3. 

c 

KZ 


I^FUT 

ROW 

DIMENSION 

OF Z 

IN CALLING 

PROGRAM. 

MIN=24. 

c 

KW1 


INPUT 

ROW 

DIMENSION 

OF W1 

IN CALLING 

PPCGR AM 

. MIN= 18 . 

c 

KW2 

=: 

INPUT 

ROW 

DIMENS ION 

OF W? 

IN CALLING 

PROGRAM 

. MIN=18. 


C 

C NE&RC-. e y .anation 

C 1 = DIMENSION SIZE EXCEEDED (KZ.KWl ,KW2) . 

C 2 = NAMEM IMPROPER LY DEFINED. 

C 

NERRCR=2 

IE (K2 .LT. 24 .OR. KW1 . LT . IE .C.E . KW2 .LT. IP) GO TO 999 
DO c J=l,?<- 
DC 5 I -1 , ?*. 

5 2 ( I , J 1 = C-.C 

IF ( NAME M .EC. 6BM] ) Cr 1C 11G 

IR (NAMEM .tO. 6HM2 ) GO TO 110 

NERR0R=2 
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GC TG 99«> 

C 

110 DO 112 1=1*3 

CWfl.l) = CJ(I,1I 
€ W < I * 1 ) = F .1(1, II 
CW (1,2) = CJ (1*2) 

EW( 1*2 I = EJ(I,2) 

CW (1*3) = C J( 1*3) 

112 EW( 1,3) = E J( I *3 ) 

CALL MAS 2 (CW , E W, TMAS *R0,NAMEM,W1 *W2, W3 *3*3 ,KW1 ,KW2 ,10) 

CALL REVADD ( .5* Wl, iVl, 1V1* 2, 16*18*24*24* KW1*K2) 

DC 113 1=1,3 
CW(I*1) - CJ(I,1) 

EW ( I » 1 ) = E J( I *1 ) 

CW ( 1,2) = CJ(I*3) 

EK (1*2) = FJ(I*3) 

CW (1,3) = CJ(I,4) 

113 EW (1,3) = E J( 1*4) 

CALL MA Z2 (CW,EW, TM AS ,R(*NAMEM,W1 ,W2,W3,3,3 ,KWl ,KW2,10) 
CALL REVADD ( . 5,W1 , 1V2, IV2, Z, 16,18,24,24* KW1,KZ) 

CO 114 1=1,? 

CW(I,1 ) = CJ(I,1 ) 

EW (1*1) = r J( I * 1 ) 

CW (1*2) = C J( 1*2) 

EW ( 1,21 = E J( 1,2) 

CW (1*3) = C J ( I *4 | 

114 EW (1*3) = E J ( I ,4 ) 

CALL MAS2 (CW,EW,TMAS,Rr*NA > ’FM*Wl *W2* W3 *3 *3 *KW1 *KW2 ,10) 

CALL REVADD ( .5,W1 , IV3, IV3,Z, 18,18*24,24, KWl,KZ) 

DO ! 1 5 1=1,3 
V (1,1) = C J( I *2 ) 

E W ( I , 1 ) = E J( I *2 ) 

C W (1*2) = C J ( I *3 > 

EW ( 1*21 = E J( I ,? ) 

CW ( 1*3) = C J ( I ,4 ) 

115 E W ( I *2 ) = EJ(I,4) 

CALL MAS? (CW,RW,TMAS,Rr,NAKEH,Wl,W2,W3,3,3,KWl,KW2,10) 
TALL kEVADD ( .5,W1,IV4,IV4*Z, ie*18,24,24, KW1,KZ) 

RETURN 

C 

999 CALL Z26CKB (5HKAS3 ,NERRC.'R) 

END 
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SUBROUTINE MAS3A I CJ ,E J. TMAS ,RO,NAME M,Z,W1 ,W2*KC J,KE J,KZ *KW1 ,KW2 I 
DIMENSION C Jf KCJ • 1) • EJlKEJ.l), Z(KZ*1), Wl(KWl,l), W21KW2,1) 

C 

C SUBROUTINE TC CALCULATE FINITE ELEMENT... 

C MASS MATRIX 

C FOR A RECTANGULAR SMEAR PANEL ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C MASS MATRIX IS IN GITRAL COORDINATE DIRECTIONS. 

C GLOBAL COORDINATE ORDER IS 
C <U,V,W) JOINT 1, THFN JOINT 2, 3, 4. 

C WHERE U,V,W APE TRANSLATIONS. 

C EULER ANGLE CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

C CALLS FORMA SUE ROUTINES M3C1,ZZE0MB. 

C DEVELOPED BY RL KC'HLEN. APRIL 1974. 

C LAST REVISION BY RL WOHLEN. HAY 1976. 

C 

C SUBROUTINE ARGUMENTS 

C CJ = INPUT MATRIX OF GLOBAL X,Y,Z COORDINATES AT PANEL JOINTS. 

C °PUS 1,2.3 CORRESPOND TO X,Y,Z COORDINATES. 

C COLS 1,2,3,* CCRRESP6ND TO JOINTS 1,2, 3, 4. SIZE (3,4). 

C EJ = INPUT MATRIX OF EULFP ANGLFS CDEGPFES) AT PANEL JOINTS- 

C ROWS 1,2,3 CORRESPOND TO GLOBAL X,Y,Z PERMUTATION. 

C CrLS 1,2 ,3,4 CORRESPOND TO JOINTS 1,2, 3, 4. SIZE(3,4). 

C TKAS = INPUT EFFECTIVE MA SS - THICKNESS . 

C RO = INPUT MASS DENSITY. 

C NAMEM = INPUT TYPE OF MASS MATRIX WANTED. 

C = Ml , LUMPED. 

C Z = OUTPUT MASS MATRIX. SIZE(12,121- 

C W1 = INPUT WORK SPACE MATRIX. SIZE(12,12I. 

C W2 = INPUT WORK SPACE MATRIX. SIZE(**,**I. 

C KCJ = INPUT ROW DIMENSION OF CJ IN CALLING PROGRAM. MIN=3. 

C KEJ “ INPUT ROW DIMENSION OF EJ IN CALLING PROGRAM. MIN'=3. 

C KZ = INPUT ROW DIMENSION OF Z IN CALLING PROGRAM. MIN=72. 

C KW1 = INPi T ROW DIMENSION OF Wl IN CALLING PROGRAM. MIN=^2. 

C KW2 = INPUT ROW DIMENSION OF W2 IN CALLING PROGRAM. MIN- *. 

C 

r NPFppp EXPLANATION 

C 1 = DIMENSION $ I i. F EXCEEDED IK Z ,KW1 ,KW? I . 

C 2 = NAMEM IMPRCPFPLY DEFINED. 

C 

NERRCP=1 

IF (KZ .LT. ] 2 .OF. KW1 .IT. 12 .OP. KW2 .LT. 0) GO TO 999 
S L 1 2 = SCRT(tCJ(l v 2)-CJIltl)l**2 ♦ IC J( ?, 2 |-C J { 2, 1 1 ) **2 

* «■ (CJ(3,2)-CJ(3,1)>**2) 

S L 1 * = SCR T ((LJ(1,4)-CJ(1, 111**2 ♦ <CJ(?,4)-CJ(2,1))**2 

* ♦ (CJ(3,4)-CJ(3,1) )**2 1 

IF ( NAM r M .^C. 6HM1 ) GO TO 110 

GC TO poo 
C 

C LUMPED. 

11C DP 112 J-1,12 
DO 1 1 T 1= 1,12 
112 Z I i » J ) = r-.u 

CALL M3C 1 (Ft!?,SLl*.,TMAS,FO,W!,KWU 
DO 1 1 B T W = 1 , A 


NERR0R=2 
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12 = 3*(IW-1) 

Z(I2*1,I2*1) = W1(IW,IW) 
2(12*2,12*2) = W1(IW,IW) 

115 2(12*3,12*3) = W1(1W,IW) 
RETURN 

999 CALL 2260MB (6HMAS3A ,NERR0R) 
END 
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SUBROUTINE PRESS (C J, T,N JN,NCOL,KC J,KW) 

DIMENSION CJ(KCJ,1),T(KW,1) 

DIMENSION A(8,6I,JNM(3,42) VN(3),C(3,9) ,IV(3),JV(9) 


*** SUBROUTINE TO CALCULATE FLUID ELEMENT PRESSURE TRANSFORMATION 
*** MORE DESCRIPTIVE CCM'.tNT CARDS TO BE ADDED AT A LATER DATE. 
*** DEVELOPED EY CA P L EODLEY, OCTOBER 1974. 


LAST 

REVISION 

BY C. S 

BCD LEY. 

NOVEMBER 

1974. 


DATA JNM 

* 1 y 2 , 3 , 

/ 

2 *4,3, 

3,4,1, 

1*4,2, 

1,2,3, 

6,5,4, 

* 

2,0,3, 

2,5,6, 

4,5,2, 

4,2,1, 

3,6,4, 

3,4,1, 

* 

3,5,6, 

3,2,5, 

4,5,1, 

1,5,2, 

1,3,6, 

1,6,4, 

* 

1*5*2* 

5,6,2, 

5,8,7, 

5,7.6, 

4,7,8, 

4,3,7, 

* 

1,2,4, 

2*3*4, 

1*4,5, 

4,8,5, 

2,6,7, 

2,7,3, 

* 

1,5,6, 

1,6,2, 

5, 8,6, 

6,8,7, 

3,7,8, 

3,8,4, 

♦ 

1,2*3, 

1,3,4, 

1*8,5, 

1 ,4,8, 

2,6,3, 

6,7,3 


CALL ZERO f T »NJN,NCOL,KW ) 

LO = 18 
NTF = 24 

IF (NJN .EQ. £) GO TO 5 
LO = 4 
NTF = 14 

IF (NJN • EC • fc) GO TO 5 
LO = 0 
NTF = 4 
5 CONTINUE 
C 

DO 20 N=1,NTF 
LOC = N + LO 
J1 = JNM(ltLCC) 

J2 = JNM(?,LOD 
J3 = JNM ( 3»LCC 1 

VN (1 )= (C J (2* J2 )— C J( 2t J1 J1MCJC3,J3)-CJ(3,J1 )) 

* “(C J (3 , J? )-CJ (3* Jl ) )*-(CJ (2,J3)~CJ(2, Jl) ) 
VN(2)=(CJ(3 t J2)-C J(3» J1 ) )*(CJ(1,J3)-CJ(1.J1 i) 

* -(CJ(1 yJ?)-CJ(ly Jl l!*(CJ(3yJ3)-CJt3y Jl M 
VN ( 3 1 - ( C J (1 y J? )— CJ( 1»J1 ) ) * ( C J (2,J3)— CJ(2,J1 ) ) 

* -CCJ(2yJ2)-CJ(2,Jl »l*(CJ'l,J31-CJ(ly Jl )1 

«C = SORT (VN(1»*VN( 1) + VN(?)*VN(2) + VN(3)*VN(3V> 
DO 25 1 = 1 ,3 
25 VN(I) = VMII/AC 
AC = AC/4R. 

IF (LOC .LT. 6) AC = ?.*AC 
DO 3C 1=1 y? 

I V ( I ) = JNM ( I , LOC ) 

DO 30 J=! ,? 

Jl * 3*1 - 3 ♦ J 
JL = (IV( 1 ) - 1 )*3 J 
30 JV(JI) = JL 
C 

DO 25 1 = 1 ,3 
DO 35 1 = 1-? 
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IL = I + ?*(L - 1) 

DC 35 J=1 »3 
F = 1. 

IF CL .EC. J) F = 2. 

35 CtJ,IL) = F*VN II) 

CALL REVADD ( AC ,C , I V* JV ,T ,3 ,9 t NJN*NCOL, 3,KW) 
20 CONTINUE 

DO 4-0 I=1*NJN 
DC 40 J=I »NJN 
A ( I r J ) = 0. 

DO 40 K=1,NC0L 

40 A(I*J) = A < I » J ) ♦ T(I ,K)*T(J,K| 

CALL INV1NP (A,A,NJN,8) 

CALL MULTB (A,T,NJN,N JN,NC0L*8»KW) 

RETURN 

END 
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C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


CROUTINE QUAD C XYZ, JODF ,EUL*NUTEL ,N J, 

NUTMX,NUTKX,NUTBX,NUTLT, NUTST, 

- W,T,S,KX,KJ,KE,KW) 

DIMENSION XY2 (KX,1 > , J DOF < K J, 1 ) ,EUL ( KE ,1 ! ,W (KW, 1 ) , T(KW, 1 ) ,S ( KW ,1 ) 
DIMENSION CJ(3,4),EJ(3,4) ,IV1 C 24) 

DATA NAHEL/fcHCUAD /* NRW*NR l T /2-'+»24/, I BLNK/6H /, KCJ/3, 

DATA NIT » NOT/ 5 ,6/ 

SUBROUTINE TC CALCULATE UN OPTION! FINITE ELEMENT ... 

MASS MATRICES AND IVECS (CN NUTMX), 

STIFFNESS MATRICES (SAME AS GLOBAL LOAD TRANSFORMATION MA T R ICES ) 
AND IVECS (CN NUTKX), 

UNIT LOAD BUCKLING MATRICES AND IVECS (ON NUTBXI, (NOT YETI 
LOCAL LOAD TRANSFORMATION MATRICES AND IVECS (ON NUTLT), (NOT YET! 
STRESS TRANSFORMATION MATRICES AND IVECS (ON NUTST), (NOT YET) 

FOR COMBINED MEMBRA NE — BE ND ING QUADRIf ATERAL PLATE ELEMENTS. 

MASS, STIFFNESS, BUCKLING MATRICES ARE IN GLOBAL COORDINATE 
D1RECTIT \!S . 

GLOBAL COORDINATE ORDER IS 

CU,V,W,P,Q,P) JOINT 1, THEN JOINT 2, 3, 4. 

WHERE U,V,W ARF TRANSLATIONS AND P t O,R APE ROTATIONS. 

IVEC GIVES ELEMENT DCF INTO GLOBAL DCF. EXAMPLES... 

IVEC<6)=834 PLACES ELEMENT DCF 6 INTO GLOBAL DOF 834. 

I VEC ( ?- ) =0 OMITS ELEMENT DOF 3 FROM GL 03 AL DOF. THIS CONSTRAINS 

ELEMENT DOF 3 TO ZERO MOTION. 

GLOBAL LOAD TRANSFORMATION MATRIX RELATES LOADS AT QUAD VERTICES IN 
GLOEAL COORDINATE DIRECTIONS TO DFFLcCTI CN S IN THE GLOBAL COORDINATE 
DIRECOICNS . 

ROW OR D F R IN GLOBAL LOAD TRANSFORMATION MATRIX IS 
(PU,PV,PW ,MP,MO,MR ) JOINT 1, THEN JOINT 2,3,4. 

WHERE P IS ■ OFCE AND M IS MOMENT. 

LOCAL LOAD TRANSFORMATION MATRICES. RELATES LOAD AT OUAD VERTICFS IN 
LOCAL COORDINATE SYSTEM TC DEFLECTIONS IN THE GLOBAL COORDINATE 
DIRECT IONS . 

STRESS TRANSFORMATION MATRICES PFLATES STRESS AT CUAD VERTICES IN 
LOCAL COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE 
DIRECTIONS . 

DATA ARRANGEMENT ON NUTMX, NUTKX, NUTBX, NUTLT , NUTST FOR EACH FINITE 
ELEMENT IS ( W-M,K,P ,LT,S T) 

WRITE (NUTWX) NAMFW ,N EL ,NR, NC ,NAME L , ( IB INK, 1= 1 , 5) , 

( ( W ( I , J ) , I = 1 , NR ) , J= 1 , NC ) , ( I VEC ( I ) , I =1 , NC ) 

CALLS FORMA SUE ROUTINES MAS 3, PAGFHD, STF3, ZZPOMB. 

DEVELOPED PY WA EFNFIFLD, CS BCDIFY, RL WOHLEN. MARCH 1973. 

LAST REVISION FY RL WOHLEN. MAY 1976. 


*****:?***Vv**:e#*:********VSS-************ ****************************** ****** 


INPUT DATA READ IN THIS SUBROUTINE FROM NUTEL. 
READ FROM CARDS . 

NAMEM »NAM TK ,NAM C LT , NA M r ST ,NAMEF 
P 0 , F , A N U 

TMASC,TMEMO,TFENC 

20 NEL,J1»J? , J3, J^,TMA5.V ,TMtMV,I6ENV 
IF ( J 1 .AC'. 0 ) P LTUEN 

GO TC 2 0 


IF NUTEL = 5, DATA IS 

4' PR MAT (5(46, 4X) 
FORMAT ( 3 ( 5X,E 10) ) 
FORMAT ( 3 ( 5 X , E 1 0 ) ) 
FOfMAT (515,. BIO 
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C DEFINITION OF INPUT VARIABLE 5. 

C NAMEM = TYPE OF MASS MATRIX WANTED. 

C = Ml * DIAGONAL LUMPED. OVERLAP AVERAGE C'F FOUR TRIANGLES. 

C = M ?» CONSISTENT. OVERLAP AVERAGE OF FOUR TRIANGLES. 

C = 6H OR 6 HNDM AS S » NO MASS MATRIX CALCULATED. 

C NAMEK = TYPE OF STIFFNESS MATRIX WANTED. 

C = Kl« OVERLAP AVERAGE OF FOUR TRIANGLES. 

C = 6H OR 6HN0STIF, NO STIFFNESS MATRIX CALCULATED. 

C NAMELT = IDENTIFICATION NAMF FOR LOAD TRANSFORMATION MATRICES. 

C = fcH OR 6HN0LCAD* NO LOAD TRANS FORMAT IONS CALCULATED. 

C NAMEST = IDENTIFICATION NAME FOR STRESS TRANSFORMATION MATRICES. 

C = fch OR 6HN0STRS, NO STPESS TRANSFORMATIONS CALCULATED. 

C NAMED = T YP F OF BUCKLING MATRIX WANTED. 

C = 6H CR 6HN0EUCK, NO BUCKLING MATRIX CALCULATED. 

C RO = MASS DFNSITY. 

C E = YOUNGS MODULUS OF ELASTICITY. 

C ANU = POISSONS PATIO. (E/2G)-1. 

C TMASC ^ EFFECTIVE MASS THICKNESS * (CONSTANT). 

C TMASV = EFFECTIVE MASS THICKNESS, (VARIASLF). 

C IF .LF. C.t TMASC IS USED. 

C TMEMC = EFFECTIVE MEMBRANE THICKNESS, (CONSTANT). 

C TMEMV = EFFECTIVE MEMBRANE THICKNESS, (VARIABLE). 

C IF .LE. 0., TMEMC IS USED. 

C TBENC = EFFECTIVE BENDING THICKNESS, (CONSTANT). 

C TBENV = EFFECTIVE BENDING THICKNESS, (VARIABLE). 

C IF .LE G., TBENC IS USED. 

C NEL = FINITE FLFMFNT NUMBER. FOR REFERENCE ONLY, NOT USED IN 

C C ALCDLATIONS. WRITTEN ON NUTMX , ETC. 

C J1 = JOINT NUMBER AT QUADRILATERAL VERTEX 1. 

C J2 = JOINT NUMBER AT QUA DR IL A TEFAL VERTEX 2. 

C J3 = JOINT NUMBER AT QUADRILATERAL VERTFX 3. 

C JA = JOINT NUMBER AT QUADRILATERAL VERTEX A. 

C 

C EXPLANATION OF INPUT T MATS . NUMBER INDICATES CARD COLUMNS USED. 

C I = INTFGEF DATA, RIGHT ADJUSTED. 

C E = DECIMAL POINT DATA, ANYWHFP E IN FIELD. EXPONENT RIGHT ADJUSTED 

C X = CARO COLUMNS SKIPPED. 

0 **************************** **** ************************************** *** 

c 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C XYZ = MATRIX OF JOINT GLOBAL X,Y,Z LOCATIONS. ROWS CORRESPOND 
C Tr JOINT MJMEFFS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 

C X , Y , Z LOCATIONS RESPECTIVELY. SIZF(NJ,3). 

C JDC-F = MATRIX CF JOI‘7 GLOBAL DEGREES OF FREEDOM. ROWS CORRESPOND. 

C TP JP1NT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 

C TRANSLATION DCFS ANT COLUMNS A,S,6 CORRESPOND 10 THE JOINT 

C ROTATION DOFS. SIZE(NJ,fc). 

C EUL » MATRIX rF JOINT EULER ANGLES (DEGREES). ROWS CORRESPOND 
C TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE 

C CLOB/L X , Y , Z PERMUTATION. SIZF(NJ,3). 

C NUTEL = LOGICAL NUMFFR OF TaF>F CONTAINING F LEM.FNT INPUT DATA FOR 
C THIS SUFROUTINE. IF NUTEL = b , DATA IS READ FROM CARDS. 

C NJ = NUMBER Of JCINTS 0 R ROWS IN MATPICES (XYZ), (JDOF), (EUL). 

t NUTMX - wCGICAL f UMF E r T UTILITY TAPE ON WHICH ELEMENT 
C MASS MaTRICFS AND IVFCS ARE OUTPUT. 
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NUTKX 


NU1BX 


NUTLT 


NUTST 


W 

T 

S 

KX 

KJ 

KE 

KW 


NUT MX MAY EE ZERO IF MASS MATRIX IS NCT FORMED. 

USES FORTRAN READ, WRITE. 

LOGICAL NUMBER CF UTILITY TAPE CN WHICH ELEMENT 
STIFFNESS MATRICES (SAME AS GLOBAL LOADS TRANSFORMATION 
MATRICES) AND IVELS ARE OUTPUT. 

NUTKX MAY EE ZERO IF STIFFNESS MATRIX IS NOT FORMED. 

USES FORTRAN READ, WRITE. 

LOGICAL NUMBER OF UTILITY TAPE CN WHICH ELEMENT UNIT LOAD 
BUCKLING MATRICES AND IVECS ARE OUTPUT » 

NUT EX MAY EE ZERO IF BUCKLING MATRIX IS NOT FORMED. 

USES FORTRAN READ, WRITE. 

LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT LOCAL 
LOAD TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. 

NUTLT MAY BF ZERO IF LOAD TRANSFORMATIONS ARE NCT FORMED. 
USES FORTRAN READ, WRITE. 

LOGICAL NUMBER CF UTILITY TAPE ON WHICH ELEMENT 
STRESS TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. 

NUTST MAY BE ZERO IF STRESS TRANSFORMATIONS ARE NOT FORMED. 
USES FORTRAN READ, WRITE. 

MATRIX WORK SPACE. MIN SIZE (24, 24). 

MATFIX WORK SPACf. MIN SIZE(24,24). 

MATRIX WORK SPACE. MIN SIZE ( 2^,24 ) . 

ROW DIMENSION OF XYZ IN CALLING PROGRAM. 

ROW DIMENSION CF JDPF IN CALLING PROGRAM. 

ROW DIMENSION OF EUL IN CALLING PROGRAM. 

ROW DIMENSION OF W, T, AND S IN CALLING PROGRAM. MIN=?4. 


NFRRCR EXPLANATION 

1 = JOINT NUMBER GREATER THAN NUMBER OF JOINTS. 

2 = MASS MATRIX FORM D , NUTMX .LF. ZERO. 

3 = STIFFNESS MATPIX FORMED, NUTKX .LE. ZERO. 

4 = LT MATRIX FORMED, NUTLT .LF. ZERO. 

5 = ST MATPIX FORMED, NUTST .LE. ZERO. 


1001 FORMAT 

1002 FORMAT 

1003 F f t MAT 

2001 FORMAL 

* 

2002 FORMAT 
* 

2003 FORMAT 

* 

* 

* 

* 

* 

* 

* 

* 

* 

2004 FO-MA': 

2005 FORMAT 


it (AG,4X) ) 

(3 If X'FIO.OM 
(5 IE ,3EK .C ) 

( / /2 5 X 40R1NPUT DATA FOR COMBINED MEMBRANE -BENDING 
2QH C.UADR IL ATER AL PLATE ELEMENTS) 

(//2LX 4CH1NPUT DATA FOR COMBINED MEMBRANE -BENDING 

41H CUADRILATFBAL PLATE ELFhFNTS (CONTINUED)) 

(/ I3X7HMASS = A6, 13X7HSTIF = A6, feX 1 3HL0AD TRANS = A6, 
3X1EFSTRESS TRANS = A6, 3X1 1H FUCK LING = A6, 

/ 15X4PF.0 = E 10.3, 13X3HE = E10.3, 

/ KX4HT (MASS ) = F10.3, 12XAHNU = E10.3, 

/ 32X13HT (MEMBRANE ) = FI0.3, 

/ 33X17HT (BENDING ) = f 1 0 . 3 , 

/ / 1 2 X 7HFLEMENT EX 7H JOINT 1 c x 7 H JO I NT 2 5X 7HJ0INT 3 
fX 7h Jf 1 NT A EX 7Hl (MASS ) 6X 11HT (MEME PANE ) 

E> I OAT (BENDING) 

/12X EB'NDMBE c <-b A 3 ( 5X 10H ( VAR I ABL F ) ) ) 

(12X 5(I5,7X),3((lL._,EX) ) 

(i ?x . n c ,7x ) ) 


c 

C READ AND WRITE FINITE c L B P t f Y DATA. 
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NLINF = 0 
CALL PAGE HD 
WRITE ( NOT ,2001 ) 

READ ( HUT EL * 1 001 ) NAMFM ,NAKEK,NAMFLT,NAMFST, NAMED 
READ ( NUT E L » ICC 2 ) RO, C,ANU 
READ (NUTFL » 1002 ) TMA SC ,TMFMC ,TPFNC 
WRITE (NOT, 2003) N ME M, NAMEK, NAME LT ,N AM EST »NAMEB» 

* RC »E » TMA SC »ANU»TMEMC,TBENC 
20 READ (NUTEL,10u3) NFL ,J 1 * J2,J3 ,J4,TMASV,TMEMV, TBENV 

NC' THIK = 1 

IF (TMASV.LE.O. .AND. TMEMV.LE.C. .AND. TBENV. LE.O.) NO fHIK=0 
IF (J1 .LE. 0) RETURN 
NLINE = NLINE + 1 
IF (NLINE .LE. 4?) DC TP 30 
CALL PAGE HO 
WPITE INC1,2CC2) 

WRITE (NOT, 200?) NAME M»NAMEK,NAMELT ,N«MFST *NAMEB, 

* PC,E ,TMASC,ANU,TMEMC,TPFNC 
NLINF =0 

30 IF (NO THIK.EQ.l ) 

*WR I TE (NCT,?004) NELt J1 ,J2,J3, J4,TMASV,TMEMV, TBENV 
IF » NO . THIK. E 0.0) WRITE (NCTt2005) NE Lt J1 » J2 , J3 » J4 

NERR0R=1 

IF (Jl.GT.VJ .PR. J2.GT.NJ .OR. J3.GT.N J .OR. J4.GT.NJ> Gu TO 999 

SET THICKNESS 1 S. 

TMAS = TMASC 
TMEM = TM EMC 
TEEN = TFENC 

IF (TMASV.GT.C.) TMAS =7 MASV 
IF (TMEMV.GT.O. ) TMEM =TMEMV 
IF (TPFNV.C-T.O.) TP EN =TEF.NV 

FORM FINITE ELEMENT CCORDINAT LOCATIONS, FULER ANGLES* REVADD iVcc. 
DC 42 I- 1 ,3 
CJ(I ,1 ) = XV? (J1 ,I> 

CJ(I?> = XV Z ( J2 * I > 

C J ( I , . ' = XYZ (J3,I) 

CJI!/ » = *YZ( 14*1) 

FJd.il = -rm i * i > 

E J ( T ■. / ) = r • 1 ) 

F J ( I *3 ) - .-',1) 

42 FJd,4) r. »<L ( J4 , 1 ) 

on ^ i = i 

IVl (I 1 = JDPF ( J1 ,1 ) 

1VI (I+fc ) = JDPF ( J2,I ) 

1 VI (1 + 12) = JDCF ( J3 *1 ) 

44 I VI (I + 1S) - JDC F ( J4 * I ) 

C 

C FORM MASS MATRIX (W) . 

IF AMS,' .K . 6H .PR. NAMFM /(*, 6HNPMASS) GO TO 110 

CA» M A ' (C J,F. JtTM AS*FL , NAMFM, W, I. ~ , KC J,KC J ,KW ,KW »KW ) 

NERRCR=? 

IF ( MU T MX ,LF. C) GP TP RP9 

WPITE (NCTFX) MANtMtNEL ,N«W ,NE - NAME L , ( IP LNK , 1 = 1, 3) , 
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( <W(I ,J )» 1=1, NRW) ,J=1,NRW) , (I VI ( I J , 1 = 1 »NkW) 

C 

C FORM STIFFNESS ^ATRIX »W), LOCAL LOAr TRANSFORMATION MATRIX (T), 

C STRFSS TRANS FOPi-ATIO:! MATRIX (S). 

110 IF ( NAMFK .EC. 6H .OR. NAME* -EC. 6HNCSTIF) GO TO 20 

CALL STF3 ( J J, E J, TM EM ,TPEN, E ,ANU, NAMEK ,NAMEST,W »T, S ,NRST , 

* KLJ » KC J» KW ,KW» KW ) 

NEnR0R=3 

IF ( NUTKX .LF. 0) GC TO 909 

WRITE (NUTKX) NAMEK ,NFL ,NRW ,NRW,NAMEL , ( IELNK, 1=1 , 5 ) * 

* ( (VM I ,J ) , 1=1 , NRW) *J=l »NRW ) * ( IV1 ( I ) ,1 = 1 »NRW) 

IF (NAMELT ..FQ. 6H .OP. NAMELT .EQ. 6HNOLOAD) GO TO 115 

NERRCK--4 

IF (NUTLT .LE. 0) CO TO 999 

WRITE ( NU f L*i ) NAMFLT* NEL ,NRL7 »NRW»NAMEL » ( TCl.NK , 1=1 ,5 ) » 

* I ( T 1 1 }.1=1,NRLT>* J=l,NRW»,(TVl(I) *I~1,NRV) 

115 IF (NAMES! .EG. 6F; .CR. NAMES! .EC . 6NNCSTR S ) GO TO 20 

NERR0P=5 

IF (NUTS! .LE. r ) GC TO 999 

WRITE (NLTST ) NAMES T, NEL ,NRST ,NPW *N AMEL , ( JBLNK , 1= 1 , 5 > , 

* ( (S (I ,J ),I= ! s NRST) ,J = i»NP.W) , ( IV 2 1 1 1 ,1=1 ,NRW) 

GC TC 20 

C 

999 CALL 2 Z POMP (6HCUAD *N ERROR ) 

E ND 
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St -.POTT! MB RECTSP I XYZ, JDOF,EUL,NUTFL ,NJ, 

* NUTMX, N"TKX,NUTLT, NUTST, 

» W, T,S,KX,KJ,KE,KW) 

L 1MENSIDN XYZtKX, l),JDOF|KJ,l ),EUL|KE,1 ) ,WIKW, 1 1 ,TIKW,1 ) ,S IKW.l ) 
r v )WION CJ13,*> ,FJ13,*) tIVl ( 1? 1 

. * i A NAMEL/6HRECTSP/, NRW,NRLT/12,8/, IBLNK/6H /, KCJ/3/ 

i.iTA N IT, KPT/ S, 6/ 


SUt .JUTINF TC* CALCI'LA'- (ON OPTION) FINITE ELEMENT ... 

S ASS MATRUES AN T VECS ( ON NUTMX), 

R.’IFFrxESS MATRICES (SAME AS GLOBAL LOAO TRANSFORMATION MATRICES! 
ADD IVFCE (t'N NUTKX), 

LOCAL LOAD TRANSFER MATICN MATRICES AND IVECS ION NUTLT), 

STRESS TRANSFORMATION MATRICES AND IVECS ION NUTST), 

FOm RECTANGULAR SHEAR panel elements. 

KA r .S, STIFFNESS MATRICES ARE IN GLOBAL COORDINATE DIRECTIONS. 

G: CEA). COORDINATE ORDER IS 

(U,V,W) JOINT 1, THEN JOINT 2, 3, A. 

WHERE U,V,W APE TRANSLATIONS. 

IVEC GIVES ELEMENT OOF INTO GLOBAL DOF. EXAMPLES... 

IV r 'IM=fc3A PLACES ELEMENT DCF 6 INTO GLOBAL DOF 83A. 

IViClD-O OMITS ELEMENT DOF 3 FROM GLOBAL DOF. THIS CONSTRAINS 
ELEMENT DOF 3 TO ZERO MOTION. 

GLOBAL LPAD TRANSFORMATION MA.TFIX RELATES LOADS AT PANEL VERTICES IN 
GLOBAL COORDINATE DIRECTIONS TC DEFLECTIONS IN THE GLOEAL COORDINATE 
DIRECTIONS . 

ROW ORDER IN GLOBAL LOAD TRANSFORMATION MATRIX IS 
|PU,PV,PW) JOINT 1, THEN JOINT 2, 3, A. 

WHERE P IS FORCE. 

LOCAL LOAD TRANSFORMATION MATRICES RELATES LOAD AT PANEL VERTICES IN 
LOCAL COOPDINhTE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE 
DIRECTIONS. 

STRESS TRANSFORMATION MATRICES RELATES PANEL SHEAR STRESS ICONSTANT) 
IN LOCAL COORDINATE SYSTFM TO DEFLECTIONS IN THE GLOBAL COORDINATE 
DIRECTIONS . 

DATA ARRANGE MB NT ON NUTMX, NUTKX , NUTLT , NUTST FOR EACH FINITE 
ELEMENT IS C W=M,K ,LT ,ST) 

WRITE (NUTWX) NAMEW ,N EL ,NR,NC ,NAMEL , ( IB LNK, 1=1 ,5! , 

( »WII,J),I=1,NR),J=1 ,NC ), I IVEC 1 1 )*I— 1»NC) 

CALLS FORMA SUBROUTINES MA S3 A, PAGEHD, STF3A, ZZEOKB . 

DEVELOPED PY P L WPHLEM. APRIL 19?A. 

LAST REVISION BY WA BENFIElD. MARCH 1°76. 


****?***&* **$* *£¥******$$*$* ** ***$43:*** **¥ A* £$$**** * 


INPUT DATA READ IN THIS SUBROUTINE FROM NUTEL. 
READ FROM CARDS. 

N AMEM ,N AMBK ,NAME LT, NAME ST 
PO,G 

TMA.S , TSTF 

20 NEL, JI > J2 ,J3,JA 

IF (Jl .Ft. 0) RETURN 
GO TO 20 


IF NUTEL = 5, DATA IS 

FORMAT 1*1 A6, AX) 
FORMAT I215X,E10> ) 
FORMAT |2|5X,E10)) 
FORMAT 1515) 


c -a 


DEFINITION OF INPUT VARIABLES. 

NAMEM = TYPE OF MASS MATRIX WANTED 
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NAMEK 


NAMELT 

NAME ST 

RO 

G 

TMAS 

TSTF 

NEL 

J1 

J2 

J3 

J4 


= Ml, DIAGONAL LUMPED. 

= M2, CONSISTENT. 

= EH no 6 HN PM ASS, NO MASS MATRIX CALCULATED. 

TYPE OF STJFr-NFSS MATPIX WANTED. 

= Kl, LINEAR DISPLACEMENT (CONSTANT STRAIN). 

= 6 H OR 6HN0STIE, NO STIFFNESS MATRIX CALCULATED. 

IDENTIFICATION NAME FOR LOAD TRANSFORMATION MATRICES. 

- 6H OR 6HNOLPAO, NO LOAD TRANSFORMATIONS CALCULATED. 

IDENTIFICATION NAME FOR STRESS TRANSFORMATION MATRICES. 

= 60! OR 6KN0STRS, NO STRESS TRANSFORMATIONS CALCULATED. 

MASS DENSITY. 

S>-F MODULUS OF ELASTICITY. 

EFFECTIVE MASS THICKNESS. 

EFFIOTTVE STIFFNESS THICKNESS. 

FINITE ELEMENT NUMBER. FOR REFERENCE ONLY, NOT USED IN 
CALCULATIONS- WRITTEN ON NUTMX, ETC. 

JOINT NUMBER AT PANEL VERTEX 1. 

JOINT NUMBER AT PANEL VERTEX 2. 

JOINT NUMBER AT PANEL VERTEX 3. 

JOINT NUMBER AT PANEL VERTEX 4. 


EXPLANATION OF INPUT FORMATS. NUMBER INDICATES CARD COLUMNS USED. 

I = INTEGER DATA, RIGHT ADJUSTED. 

E = DECIMAL POINT DATA, ANYWHERE IN FIELD. EXPONENT RIGHT ADJUSTED. 
X = CARD COLUMNS SKIPPED. 


SUBROUTINE ARGUMENTS (ALL INPUT) 

XYZ = MATRIX OF JOINT GLOBAL X,Y,Z LOCATIONS. ROWS CORRESPOND 
TC JOINT NUMBFRS. COLUMNS 1,2.3 CORRESPOND TO THE JOINT 
X , Y , 2 LOCATIONS RESPECTIVELY. SIZE(NJ,3). 

JDCF = MATRIX OF JOINT GLCEAL DEGREE 1 OF FREEDOM. ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 
TRANSLATION DOFS AND COLUMNS 4,5,6 CORRESPOND TO THE JOINT 
ROTATION DOFS. SI2F(NJ,6l. 

EUL - MATPIX OF JOINT EULEP ANGLES (DEGREES). ROWS CORRESPOND 
TO JOINT NUMPFRS. COLUMNS 1,2,3 CORRESPOND TO THE 
GLOBAL X ,Y ,Z PERMUTATION. S1ZF(NJ,3). 

NUTEL = LOGICAL NUMBER OF TAPE CONTAINING ELEMENT INPUT DATA FOR 
THIS SUBROUTINE. IF NUTEL = 5, DATA IS READ FROM CARDS. 

NJ = NUMBER OF JOINTS OR ROWS IN MATRICES (XYZ), (JDOF), (EUL). 

NUTMX = LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 
MASS MATRICES AND IVECS ARE OUTPUT. 

NUTMX MAY BE ZERO IF MASS MATRiX IS NOT FORMED. 

USES FORTRAN PE AD, WRITE. 

NUTKX = LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 

STIFBNFSS MATRICES (SAME AS GLOBAL LOADS TRANSFORMATION 
MATP1CFS) AND iVECS APE OUTPUT. 

NUTKX MAY r 'F ZERO IB STIFFNESS MATRIX IS NOT FORMED • 

USES FORTRAN °E AD , WRITE. 

NUTLT = LOGICAL NUMBfP OF UTILITY TAPE ON WHICH ELEMENT LOCAL 
LOAD TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. 

NUTLT MAY BF 7EFO IB LOAD TRANSFORMATIONS ARE NOT FORMED. 
USES FORTRAN PFAD, WRITE. 

NUTST = LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 



on on 


RECTSP — 3/ 4 


C STRESS TP AN S F OR MA TI ON MATRICFS AND IVECS ARE OUTPUT. 

C NUTST MAY PE ZERO IF STRESS TRANSFORMATIONS ARE NOT FORMED. 

C USFS FORTRAN READ, WRITE. 

C W = MATRIX WORK SPACE. MIN SIZF(12»12). 

C T = M AT r IX WCFK SPACE. MIN SIZEI12,12). 

C S = MATRIX WORK SPACE. MIN S1ZE(12,12). 

C KX =• ROW DIMENSION OF XYZ IN CALLING PROGRAM. 

C KJ = ROW DIMENSION OF JDOF IN CALLING PROGRAM. 

C KE = f PW DIMENSION OF EUL IN CALLING PROGRAM. 

C KW = ROW DIMENSION OF W, T, AND S IN CALLING PROGRAM. MIN=12. 

C 

C N ERROR EXPLANATION 

C 1 = INPUT JOINT NUMPEP EXCEEDS MAXIMUM ALLOWABLE NUMBER OF JOINTS. 

C 2 = NUTMX NON POSITIVE. 

C 3 = NUTKX NON POSITIVE. 

C 4 = NUTLT NON POSITIVE. 

C 5 = NUTST NON POSITIVE. 

C 

1001 FORMAT (4(A6,4X)) 

1002 FORMAT { 2 (5X, E10 .0 ) ) 

1003 FPP MAT €513) 

2001 FORMAT <//3£X 47H INPUT DATA FOP RECTANGULAR SHEAR PANEL ELEMENTS) 

2002 FORMAT (//32X ^7H INPUT DATA FOR RECTANGULAR SHEAR PANEL ELEMENTS 

* 12H (CONTINUED)) 

2003 FORMAT (/ 14X7HM ASS = A6, 1 4X7HSTIF = A6, 11X13HL0AD TRANS = A6 , 

* 8X15HSTRESS TRANS = A6, 

* / 16X4HRD = E 10.3 , 14X3HG = E10.3, 

* / I1XPHT (MASS) = E10.3, 8X9HT(STIF) = E10.3, 

* //1PX7HELEMENT 13X7HJOINT 1 13X7HJ0INT 2 13X7HJ0INT 3 

* 13X7H JOINT 4 / 18X6HNUMBE R ) 

200* FORMAT (1FX,S (15,lbX» ) 

READ AND WRITE FINITE ELEMENT DATA. 

NLINE = 0 
CALL PAGE ED 
WRITE (NOT, 2001) 

READ ( NUT FL, 10C1 ) NAM EM ,N AMEK ,NAMELT,NAMEST 
READ ( NUT FL, 1002 ) RO*G 
READ (NUT EL , 1 002 ) TMAS»TSTF 

WRITE ( NOT, 2003) NAMEM,NAMEK, NAMELT ,NAMEST,RO,G,T MAS.TSTF 
20 READ ( NUT FL, 1003) MEL ,J1»J2»J3»J4 
IF (J! .LF. 0) RFTURN 
NL INF = NLINE + 1 
IF (NLINF .LF. 42) GO TO 30 
CALL PAGFED 
WRITF (NOT, 2002) 

WRITF (NOT, 2003) NAME M,NAMFK, NAMELT, NAMES!, RO,G,TMAS,TSTF 
NLINE = 0 

30 WRITE (NCT, 2004) NEL»JI»J2,J3»J4 

NERR0R=1 

IF (Jl.GT.NJ .OR. J2.GT.NJ .OR. J3.GT.NJ .OR. J4.GT.NJ) GO TO 999 

FORM FINITE ELEMENT COORDINATE LOCATIONS, EULER ANGLES, REVADD IVEC. 
DO *2 1-3 ,3 

cj(i,n * xyz ( ji , i ) 
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C J ( I » 2 ) = XYZ ( J2 » 1 ) 

CJ(I,3) = X Y2 ( J3 » i ) 

CJ( I ,4 ) = XYZ ( J4, I ) 

E J ( I »1 ) - FUL(J1 ,1) 

EJ(I,2) = FUL(J2,n 
E J ( I »3 ) = EUL(J3,I) 

42 E J ( I »4 ) = FUL(J4,1) 

Dr 44 1=1 t 3 

1V1 (I) = JDPF ( J 1 t 1 ) 

IV1 (1+3 ) - JDOF ( J2» I ) 

IV1 (1+6 ) = JDCF ( J?* I ) 

44 I VI (1+9 1 = JDPF(J4,1) 

FORM MASS MATRIX (W). 

IF ( NAME M .EC . fcH . OF. . NAMEM .EC. 6HNCMASSI GO TO 110 

CALL MASS A ( CJ,EJ*TMAS,RO,NAMEM,W,T,S,KGJ»KCJ*KW > K.«,KW 1 

NERR0R=2 

IF INUTMX .LE. 0) GO TO 999 

WRITE (NUTMX) NAMEM ,NFL ,NRW ,NRW,NAMEL , ( I6LNK , 1=1 , 5 ) . 

* UW(1»J ),I=1,NRW) ,J=1 ,NRW) t (IV1 (I), I=1,NRW) 

C 

C FORM STIFFNESS MATRIX (W ) * LOCAL LOAD TRANSFORMATION MATRIX (Tl* 

C STRESS TRANSFORMATION MATRIX (S). 

110 IF ( NAMFK .FQ. 6H .CP. NAMEK .EO. 6HN0STIF) GO TO 20 

CALL STF3A { CJ, EJ, TSTF ,G .NAMFK, NAMFST, W,T, S ,NRST , 

* KCJ,KCJ*KW*KW,KW) 

NERR0R=3 

IF (MUTKX .LF. 0) GO TO 999 

WRITE (MUTKX) NAMEK. ,N EL ,NRK ,NRW,NAMEL * ( IBLNK, 1 = 1, 5 ) , 

* ( (WII'J ),I=1'NRW) »J=1 *NRW), (I VI (I)*I=1,NPW> 

IF (NAMELT .EC. 6h .OR. NAMELT .EQ. 6HN0L0AD) GO TO 115 

NERR0R=4 

IF (NUTLT .LE. C) C° TO 9«9 

WRITF (NUTL1) NAMEL T, NF L* NRLT »NkW,N AMEL , ( IRLNK , 1= 1 ,5 ) » 

* ( ( T ( 1 » J ) » 1=1 t NRLT) t J= 1 »NRW) * ( IV1 ( I ) »I=1*NRW) 

115 IF (NAMES! .EC. 6H .OP. NAMEST .EQ. 6HN0STP S) GO TO 20 

NERR0R=5 

IF (NUT ST .LF. 0) GO TO 9«9 

WRITE (NO 1ST) NAMEST»NEL,NRST,NPW,NAMELtnP»-NK,I=l,5), 

* ((S(I»J)»I= I* NRST) * J = 1»NRW) *(IV1(I) ,I=1*NRW) 

GO TO 20 

C 

999 CALL ZZBOMF ( 6HR ECT SP ,NFRROR ) 

END 
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SUBROUTINE STF1A { C J ,E J, Al ,A 2, E ,N AMEK,NAME ST , S,TL,TS ,NRST» 

* KC J,KEJ,KS,KTL,KTS) 

DIMENSION Cj(KCJ, 1), EJ(KEJ,1), S(KS,li, TL(KTL,1), TSfKTS,!) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C STIFFNESS MATRIX ( S AMF AS GLOBAL LOAD TUANS FORMATION MATRIX), 

C LOCAL LOAD TRANSFORMATION MATRIX, 

C STRESS TRANSFORMATION MATRIX, 

C Fnp an AXIAL ROD ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C ROD MAY BE LINEARLY TAPERED CR UNIFORM. 

C STIFFNESS MATRIX IS IN GLOBAL COCRDINATF DIRECTIONS. 

C GLOBAL COORDINATE ORDER IS 
C (U,V,W) JOINT 1, THEN JOINT 2. 

C WHERE U,V,W ARE TRANSLATIONS. 

C GLOBAL LOAD TFaNS’OPMATI ON MATRIX RELATES LOADS AT ROD ENDS IN GLOBAL 
C COORDINATE DIRECTIONS TO DEFLECTIONS IN THE GLOBAL COORDINATE 
C DIRECTIONS. 

C ROW ORDER IN GLC'EAL LOAD TRANS FORMATION MATRIX IS 
C (PU,PV,RVT) JOINT 1, THEN JOINT 2. 

C WHERE P IS FORCF. 

C LOCAL LOAD TRANSFORMATION MATRIX PELATES LOADS AT ROD ENDS IN LOCAL 
C COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE DIRECTIONS. 
C ROW ORDER IN LOCAL LOAD TRANSFORMATION MATRIX IS 
C PXI ,RX2 

C WHERE PX IS AXIAL FORCE. 

C PXl(-), PX2 ( ♦) IS TFNSION. PXK + ), PX21-) IS COMPRESSION. 

C STRESS TRANSFORMATION MATRIX PELATES STRESS AT POD ENDS IN LOCAL 
C COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOEAL COORDINATE DIRECTIONS. 
C ROW ORDER IN STRESS TRANSFORMATION MATRIX IS 
C SIGMA— XT , SIGMA— X2 

C WHERE SIGMA IS NORMAL STRESS. 

C SXl(-), $X2(+) IS TENSION. SX1( + ), SX2(-) IS COMPRESSION. 

C EULER ANGLF CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

C CALLS FORMA SUBROUTINES ATXBAl, DCOSlA, K1A1, MULTA, ZZ50MB. 

C DEVELOPED FY RL WOHLEN. SEPTEMBER 1*572. 

C LAST REVISION BY WA BENFIELD. MARCH 1*576. 

C 

C SUBROUTINE ARGUMENTS 

C CJ = INPUT MATRIX OF GLCEAL X,Y,Z COORDINATES AT ROD JOINTS. 

C ROWS 1,2,3 CORRESPOND TO X,Y,Z COORDINATES. 

C COLS 1,2 CLRRE SPOND TO JOINTS 1,2. SIZE<3,2). 

C EJ = INPUT MATRIX OF EULEP ANGLES (DEGREES) AT POD JOINTS. 

C ROWS 1,2,3 CORRESPOND TO GLOBAL X,Y,Z PERMUTATION. 

C COLS 1,2 CORRESPOND TO JOINTS 1,2. SIZE(3,2). 

C Al - INPUT CROSS-SECTION AREA AT POD END 1. 

C A2 = INPUT CROSS-SECTION AREA AT ROD END 2. 

C E = INPUT YOUNGS MODULUS OF ELASTICITY. 

C NAMFK = INPUT TYPE OF ST IF MATRIX WANTED. 

C = Kl, CONSTANT AXIAL FORCE ASSUMED. 

C NAMEST = INPUT OPTION FOR STRESS TRANSFORMATION. 

C = 6H OP 6HN0STRS ,NO STRESS TRANS CALCULATED. 

C S = OUTPUT STIFFNESS MATRIX (SAME AS GLOBAL LOAD TRANSFORMATION 

C MATRIX). SIZE(6,6). 

C TL = OUTPUT LOCAL LOAD TRANSFORMATION MATRIX. SIZE (2,6). 

C TS = OUTPUT STRESS TRANSFORMATION MATRIX. SI ZE(NRST ,6) • 
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c 

NRST 

=: 

OUTPUT 

NUMPER OF ROWS 

IN STRESS TRANSFORMATION MATRIX. 

c 

KCJ 


INPUT 

ROW 

DIMENSION 

OF 

CJ 

IN 

CALLING 

PROGRAM. M1N=3. 

c 

KEJ 

- 

INPUT 

ROW 

DIMENS ION 

OF 

EJ 

IN 

CALLING 

PROGRAM. MIN=3. 

c 

KS 

— 

INPUT 

ROW 

DIMENSION 

OF 

S 

IN l 

CALLING 

PROGRAM. MIN=6. 

c 

KTL 


INPUT 

ROW 

DIMENSION 

OF 

TL 

IN 

CALLING 

PROGRAM. MIN=2. 

c 

C 

KTS 

= 

INPUT 

ROW 

DIMENSION 

OF 

TS 

IN 

CALLING 

PROGRAM. MIN=NRST • 


C NERRCR EXPLANATION 

C 1 = SIZE LIMITATION EXCEEDED. 

C 2 = NAMEK IMPROPERLY DEFINED. 

C 

NRST = 2 

NERROR= 

IF ( K S .LT. 6 .OR. KTL .LT. 2 .OR. KTS .LT. NR ST) GO TO 999 
PL = SCPT<tCjn,2)-CJ(l,l)l**2 ♦ <CJ(2,?)-CJ(2,1H**2 
* + <CJC3.2)-CJ(3,1) )**2) 

IF (NAMEK .f Q. 6HK1 ) GO TO 110 

NERROR= 

GO TO 9°9 

110 CALL K1A1 ( Al,A2,RL,E,TL,TS,KTL,KTS) 

C 

CALL DC OS 1.A ( C J ? E J , S» KC J, KE J,KS ) 

CALL MULTA (TL,S, 2,2,6 t KTL,KS) 

IF (NAMFST .EO. 6F .OR. NAMES! .EO. 6HN0STRS) GO TO 210 

CALL MULTA < TS »S ,NPST,2,6, KTS,KS> 

210 CALL ATXBA1 (S,TL, ?,6, KS,KTL) 

RETURN 

C 

999 CALL ZZBOMB ( 6HSTF1 A »NER POP ) 

END 


2 

TL=K 

S=DC 
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SUBROUTINE STF1B (CJ ,E J,KODE ,A1 ,A2 ,TJ1 ,T J2 ,B IZI , BIZ2 ,B IY1 ,bl Y2 , 

* P.l ,R2,CY1,CY2,CZ1,CZ2,SF,E,G,NAMEK,NAMEST, 

* S,TL,TS,NRST,KCJ,KEJ,KS,KTL,KTS) 

DIMENSION C.J(KCJ,1 >,EJ(KEJ,1 ) ,KODE 1 1 > ,S (KS , 1 ) ,7 LI KTL, 1 ) ,TS CK7 S, 1 ) 

SUBROUTINE 70 CALCULATE FINITE ELEMENT... 

STIFFNESS. MATRIX (SAME AS GLOBAL LOAD TRANSFORMATION MATRIX), 

LOCAL LOAD TRANSFORMATION MATRIX, 

STRESS TP ANSFCPMATION MATP1X, 

FOR A COMBINED AX 1AL-TCR S10N-BENDING BAR ELEMENT WITH UNRESTRAINED 
BOUNDARIES. 

BAR MAY BE LINEARLY TAPERED OR UNIFORM. 

STIFFNESS MATRIX IS IN GLOBAL COORDINATE DIRECTIONS. 

GLOBAL COORDINATE ORDER IS 

<U,V,W,P ,0,R I JOIN! 1, THEN JOINT 2 
WHERF U,V,W ARE TRANSLATIONS AND P,Q,R ARE ROTATIONS. 

GLOBAL LOAD TRANSFORMATION MATRIX RELATES LOADS AT BAR ENDS IN GLOBAL 
COORDINATE DIRECTIONS TO DEFLECTIONS IN THE GLOBAL COORDINATE 
DIRECTIONS. 

ROW ORDFP IN GLOEAL LOAD TRANSFORMATION MATRIX IS 
(PU,PV,PW,MP,MQ,MR) JOINT 1, THEN JOINT 2. 

WHERE P IS FORCE AND MIS MOMENT. 

LOCAL LOAD TRANSFORMATION MATRIX RELATES LOADS AT BAR ENDS IN LOCAL 
COORDINATE SYSTFM TO DEFLECTIONS IN THE GLOBAL COORDINATE DIRECTIONS. 
ROW ORDER. IN LOCAL LOAD TRANSFORMATION MATRIX IS 
PX1,PX2,MX1,MX2,PY1,PY2,MZ1,MZ2,PZ1,PZ2,MY1 ,MY2 
WHERE P IS FORCE AND M IS MOMFNT • 

STRESS TRANSFORMATION MATRIX RELATFS STRESS AT BAR ENDS IN LOCAL 
COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE DIRECTIONS. 
ROW ORDER IN STRESS TRANSFORMATION MATRIX IS 
PX1/AI ,PX2/A2, MX1*R 1/TJ1 ,MX2*R2/TJ2 , 
PY1/A1,PY2/A2,MZ1'*CY1/B1Z1,MZ2*CY2/BIZ2, 

PZ 1 /Al, PZ2/A2, MY 1*C21/B I Yl, MY 2*C22/RIY2 
WHERE P IS FORCE AND M IS MOMENT. 

EULER ANGLE CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

CALLS FORMA SUBROUTINES ATXBA1 ,DC0S1B,K1A1,K1B1,K 1C1 ,MULTA,ZZBOMB. 
DEVELOPFD BY RL WOHLEN. FEBRUARY 1973. 

LAST REVISION BY RL WOHLEN. APRIL 19?6. 

SUBROUTINE ARGUMENTS 

CJ = INPUT MATRIX OF GLOBAL X,Y,Z COORDINATES AT BAR JOINTS. 

ROWS 1,2,3 CORRESPOND TO X,Y,Z COORDINATES. 

COLS 1,2 CORRESPOND TO JOINTS 1,2. COL 3 CORRESPONDS 
TO RFFERFNCE POINT TO DFFIMF LOCAL XY PLANE. SIZE(3,3). 
FJ = INPUT MA1P1X OF FULER ANGLES (DEGREES) AT BAR JOINTS. 

ROWS 1,2,3 CORRESPOND TO GLOBAL X,Y,Z PERMUTATION. 

COLS 1,2 CORRESPOND TO JOINTS 1,2. SIZE (3,2). 

KODF - INPUT OPTION CODF FOP AXIAL, TORSION, BONDING 2, BENDING Y 

LOCAL STIFFNESS. IF BLANK, ALL FOUR ARE CALCULATED. 

SIZE (A). 

KODE ( 1 )= A , LOCAL STIFFNESS MATRIX IS CALCULATED 
FOR AXIAL (ALONG LOCAL X-AXIS). 

KODE (2 )= T , LOCAL STIFFNESS MATRIX IS CALCULATED 
FOR TORSION (ABOUT LOCAL X-AXIS). 

KODE (3 )=BZ, LOCAL STIFFNESS MATRIX IS CALCULATED 
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A1 

A2 = 

TJ1 

TJ2 

BIZI 

B1Z2 

B1Y1 

eiY2 

R1 

R2 = 

CY1 

CY2 

CZ1 

CZ2 

SF 


E 

G 

NAMEK = 


NAME ST = 

s 

TL 

TS 

NR ST 

KCJ 

KEJ 

KS 

KTL 

KTS - 


FCR BENDING (ABOUT LOCAL Z-AX1S). 

KC>DE K) = 6Y, LOCAL STIFFNESS MATRIX IS CALCULATED 
FOP PENDING (ABOUT LOCAL Y-AXIS). 

INPUT CROSS-SECTION AREA AT BAR END 1. 

INPUT SAME AS A1 AT PAR END 2. 

INPUT CROSS-SECTION SAINT VENANTS TORSION CONSTANT (J) IN 
JG AT PAR END 1. 

INPUT SAME AS TJ1 AT BAR END 2. 

INPUT CPC'SS-SE CT ION AREA MOMENT OF TNERTIA ABOUT LOCAL 
Z-AXIS (FOR BENDING) AT BAR END 1. 

INPUT SAME AS BIZI AT EAR END 2. 

INPUT CROSS-SECTION AREA MOMFNT OF INERTIA ABOUT LOCAL 
Y-AXIS (EOF BENDING) AT PAR END 1. 

INPUT SAME AS EIY1 AT EAR END 2. 

INPUT DISTANCE FROM LOCAL X-AXIS TO OUTER FIBER FOR 
TORSION STRESS CALCULATION AT BAR END 1. 

INPUT SAME AS R1 AT PAR END 2. 

INPUT DISTANCE FROM XZ PLANE TO OUTER FIBER FOR BENDING 

STRESS CALCULATION AT BAR END 1. LOCAL Y DIRECTION. 
INPUT SAME AS CY1 AT PAR END 2. 

INPUT DISTANCE FROM XY PLANE TO OUTER FIBER FOR BENDING 

STRESS CALCULATION AT BAR END 1. LOCAL Z DIRECTION. 
INPUT SAME AS CZ1 AT BAR END 2. 

INPUT SHAPF FACTOR ( K ) FOR SHEAR IN KAG. 

US F SF=0 .0 FOR NO SHEAR DEFORMATION IN BENDING. 

SF=1 .0 FCR A SOLID CIRCULAR CYLINDER. 

S F= . 5 FOP A THIN WALLED CIRCULAR CYLINDER. 

INPUT YOUNGS MODULUS OF ELASTICITY. 

INPUT SHEAR MODULUS OF ELASTICITY. 

INPUT TYPE OF STIF MATP1X WANTED. 

= Kl, USES K 1 A 1 FOR AXIAL, K1C1 FOR TORSION, 

K) 61 FOP BENDING. 

INPUT OPTION FOR STRESS TRANSFORMATION. 

= 6H OR 6HNCSTRS ,NC STRESS TRANS CALCULATED. 

OUTPUT STIFFNESS MATRIX (SAME AS GLOBAL LOAD TRANSFORMATION 
MATRIX). SIZE! 12,12) • 

OUTPUT LOCAL LOAD TRANSFORMATION MATRIX. SIZE(12,12). 

OUTPUT SIRFSS TRANSFORMATION MATRIX. SIZE (NR.ST»12 ) • 

OUTPUT NUMBER OF ROWS IN STRESS TRANSFORMATION MATRIX. 

INPUT PDW DIMENSION OF CJ IN CALLING PROGRAM. MIN=3. 

INPUT ROW DIMENSION OF EJ IN CALLING PROGRAM. MIN=3. 

INPUT ROW DIMENSION OF S IN CALLING PROGRAM. MIN=12. 

INPUT ROW DIMENSION OF TL IN CALLING PROGRAM. MIN=12. 

INPUT ROW DIMENSION OF TS IN CALLING PROGRAM. MIN=NRST. 


NtPPP t EXPLANATION 

1 = S1ZF LIMITATION FXCEFDED. 

2 = NAMEK IMPROPERLY DEFINED. 


NR ST = 12 

IF (K5 .LT. 12 .OR. KTL 
DO 5 J = 1 » 1 2 
DO 5 1=1,12 
TL ( I , J ) = 0.0 


NERROR= 1 

.LT. 12 .OP. KTS .LT. NRST) GO TO 999 
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5 TS(I,J) = 0.0 

RL = SORT ( ( C J ( 1*2')— CJ ( 1 * 1 ) ) + ( C J < 2* 2 )-C J ( 2 , 1) )**2 

* + ( C J( 3* 2 )-C J ( 3 » 1) )**2) 

KODFA = 1 
KODET = 1 
KODEBZ = 1 
KCIDEBY = ) 

IF (KODEU).FC.IH .AND. KC'DE ( 2 ) .E 0 .1H .AND. 

* KODt ( ?) .E0.2H .AND. KCDE (4 ) .EQ .2H ) GO TO 10 

IF (KC’DF(l) .NE. 1HA ) KODEA = 0 

IF ( KOD E ( 2 ) .NE. 1HT ) KODET = 0 

C LAST HALF OF NEXT TWO CARDS ALLOW FOR OLD DATA. INSERTED APRIL 1976. 

IF ( KPDF ( 3 ) .NE. 2HB2 .AND. K0DE(3) .NE. 2HBY ) KODEBZ = 0 

IF ( KODF ( 4 ) .NE. 2HBY .AND. K0DE(4) .NE. 2HBZ) KODEBY = 0 

10 IF (NAMEK .EQ. 6HKI ) GO TO 110 

NER R0R=2 

GO TO 999 

C 

C AXIAL = K1A1 (CONSTANT FORCE), TORSION = K1C1 (CONSTANT TORQUE), 

C BENDING - K1B1 (CONSTANT SHEAR, LINEAR BENDING MOMENT). 

110 IF (KODEA .EC. 1) CALL KlAl (A1,A2,RL»E,TL ,TS »KTL,KTS ) 

IF (KODET .ECU 1) CALL K1C1 (T J1 , TJ2 ,R 1 ,R2 ,R L,G, TL (3 ,3 ) ,TS (3 ,3 ) , 

* KTL ,KTS ) 

IF (KODEBZ .EC. 1) CALL K1B1 ( BIZI , B I Z2 ,CY1 ,C Y2 ,Al , A2 ,SF , RL, E ,G , 

* TL( 5 ,5) ,TS(5,5) ,KTL ,KTS ) 

DC 115 J=7 *8 

DO 115 1=5,6 
TL ( I ,.l ) =— TL ( I , J ) 

T S ( I , J ) = — T S ( 1 , J ) 

TL ( J , I ) =— TL ( J , 1 ) 

115 TS ( J , 1 ) =-TS(J,I) 

IF (KCDFBY .EO. 1) CALL K1P1 ( BI Y1 , PI Y2 ,C Z1 ,CZ2 , A 1 , A2 ,SF ,RL,E ,G , 

* TL(9,9) ,TS(9,9) ,KTL,KTS) TL=K 

C 

CALL DC OS IB ( C J , E J , S, KC J,KE J , KS ) S=DC 

CALL MULT A (TL,S, 12,12,1 2, KTL ,KS ) 

IF (NAMEST .FC. t>h .OR. NAMEST .EC. 6HN0STR S) GO TO 210 

CALL MULTA ( TS,S,NR$T, 12 ,12, KT<,K$) 

210 CALL ATXEA1 ( S ,Tt , 12?12, KS.kTL) 

RETURN 


C 


999 CALL ZZEOME ( 6HSTF1B ,N ERROR) 
END 
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SUBROUTINE STF? (C J, F J ,TMEM,T' FN , E , ANU ,NAMEK ,NAM EST,S ,TL ,TS,NR ST , 

* KCJ ,KEJ ,KS,KTL,KTS) 

DIMENSION CJ(KCJ,1), EJ(KEJ,1), S(KS,1), TL(KTL,1), TS(KTS*1) 

subroutine to calculate finite element... 

STIFFNFSS MATRIX (SAME AS GLOBAL LOAD TRANSFORMATION MATRIX)* 

LOCAL LOAD TRANSFORMATION MATRIX, 

STRESS TRANSFORMATION MATRIX, 

FOR A COMBINED MEME R ANE- BEND 1NG TRIANGLE PLATE ELEMENT WITH 
UNRESTRAINED BOUNDARIES. 

STIFFNESS MATRIX IS IN GLOBAL COORDINATE DIRECTIONS. 

GLOBAL COORDINATE ORDER IS 

(U,V,W,P,Q,R ) JOINT 1, THEN JOINT 2, 3. 

WHERE U,V,W APE TRANSLATIONS AND P,0,R APE POTATIONS. 

GLOBAL LOAD TRANSFORMATION MAT FIX RELATFS LOADS AT TRNGL VERTICES IN 
GLOBAL COORDINATE DIRECTIONS TO DEFLECTIONS IN THE GLOBAL COORDINATE 
DIRECTIONS. 

ROW ORDER IN GLOBAL LOAD TRANSFORMATION MATRIX IS 
(PU,PV,PW ,MP,MO,MR) JOINT 1, THEN JOINT 2,3. 

WHERE P IS FORCE AND M IS MOMENT. 

LOCAL LOAD TRANSFORMATION MATRIX RELATES LOADS AT TRNGL VERTICES IN 
LOCAL COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE 
DIRECTIONS. 

ROW ORDFR IN LOCAL LOAD TRANSFORMATION MATRIX IS 
( PX, PY ,MZ ) JOINT 1 THEN 2,3, NEXT 
{ PZ , MX *MY ) JOINT 1 THEN 2,3. 

WHERE P IS FORCE AND M IS MOMENT. 

STRFSS TRANSFORMATION MATRIX RELATES STRFSS AT TRNGL VERTICES IN LOCAL 
COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE DIRECTIONS. 
ROW ORDFR IN STPELS TRANSFORMATION MATRIX IS 

( SIGMA— X , SIGMA— Y , TAU— XY ) FOR (Z=T6EN/2) AT JOINT 1, 

THEN JOINT 2,3. 

(SIGMA-X,SIGMA-Y,TAU-XY) FOR (Z=-TBEN/2) AT JOINT 1, 

THFN JOINT 2,3. 

WHERF SIGMA IS NORMAL STRFSS AND TAU IS SH FAR STRESS. 

EULFR ANGLE CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

CALLS FORMA SUBROUTINES ATXBA1 ,DC0S2,K2A1 , K2B1,MULTA ,ZZBOMB. 

DEVELOPED BY WA PENEIELD. FEBRUARY 1973. 

LAST REVISION BY WA BENE IE LD . MARCH 1976. 

SUBROUTINE ARGUMENTS 

CJ - INPUT MATRIX OF GLOBAL X,Y,Z COORDINATES AT TRIANGLE JOINTS. 

ROWS 1,2,3 CORRESPOND TO X,Y,Z COORDINATES. 

COLS 1,2,3 CORRESPOND TO JOINTS 1,2,3. SIZE(3,3). 

EJ = INPUT MATRIX OF EULER ANGLES (DEGREES) AT TRIANGLE JOINTS. 

POWS 1,2,3 CORRESPOND TO GLOBAL X,Y,Z PERMUTATION. 

COLS 1,2,3 CORRESPOND TO JOINTS 1,2,3. SIZE(3,3). 

TMEM = INPUT EFFECTIVE MEMBRANE THICKNESS. 

TbEN = INPUT EFFECTIVE BENDING THICKNESS. 

E = INPUT YOUNGS MODULUS OF ELASTICITY. 

ANU = INPUT POISSONS RATIO. (E/2G<-1. 

NAMFK = INPUT TYPE OF ST IF MATRIX WANTED. 

= K1 , USES K2A1 FOR MEMBRANE, K2B1 FOR BENDING. 

NAMEST = INPUT OPTION FOR STRESS TRANSFORMATION. 

= 6H OP 6HN0STRS ,NO STRESS TRANS CALCULATED. 
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C S = OUTPUT STIFFNFSS MATRIX (SAME AS GLOBAL LOAD TRANSFORMATION 

C MATRIX). SIZE! 18 « 18) • 

C TL = OUTPUT LOCAL LOAD TRANSFORMATION MATRIX. SIZE (18, 18). 

C TS = OUTPUT STRESS TRANSFORMATION MATRIX. SI ZE (NRST, 18 ) . 

C NRST = OUTPUT NUMBER OF ROWS IN STRESS TRANSFORMATION MATRIX. 

C KCJ = INPUT ROW DIMENSION OF CJ IN CALLING PROGRAM. 

C KEJ = INPUT ROW DIMENSION OF EJ IN CALI ING PROGRAM. 

C KS = INPUT ROW DIMENSION OF S IN CALLING PROGRAM. MIN=18. 

C KTL = INPUT ROW DIMENSION OF TL IN CALLING PROGRAM. MIN=18. 

C KTS = INPUT ROW DIMENSION OF TS IN CALLING PROGRAM. MIN-NRST • 

C 

C NERROR EXPLANATION 

C 1 = SIZE LIMITATION EXCEEDED. 

C 2 = NaMEK IMPROPERLY DEFINED. 

C 

NRST = 18 

NE_RROR=l 

IF (KS .LT. i8 .OR. KTL .LT. 18 .OR. VTS .LT. NRST) GO TO 999 
DO 5 J=1 ,18 
DO 5 1=1,18 
TL ( I , J ) = 0.0 
5 TS ( 1 , J ) = 0.0 

SL12 = SQRT ( ( C J( 1 ,2 )—CJ (1,1)) **2 + (C J( 2,2)-CJ ( 2, 1 ) )**2 

* + (CJ(3,?)-CJ(3,1) )**2) 

SL23 = SORT ( (CJl 1 ,3 )-CJ (1 ,2 ) )**2 + <CJ( 2,3)-C J( 2,2) )**2 

* + (CJ( 3,3)-CJ(3,2) )**2) 

SL13 = SQR'I((CJ<1,3)-CJ(1,1))**2 + (CJ( 2,3)— CJ(2»1) )**? 

* + (CJ(3,3)— CJ(3,1) )**2 ) 

X3 = (SL13**2+SL12**2-SL23**2)/(2.0*SL12) 

Y3 - SQRT ( S L 1 3**2 -X 3**2 ) 

IF (NAMEK . EQ. 6HK1 ) GO TO 110 

NERR0R=2 

GO TO 9°9 
C 

C MEMBRANE = K?Al ( BOD LE Y, BENF IE LD ) , FENDING = K2B1 (BODLEY). 

110 CALL K2A1 ( SL12 ,X 3. Y3 ,TMEM, F , ANU, TL ,TS, S,KTL ,KTS,KS ) TL=K 

CALL K2B1 (SL12»X3»Y3»TBEN,E» ANU ,TL(10,10)*TS(1,10),S, 

* KTL, KTS, KS) 

DO 111 1=1,9 

II = 1+9 
DC 111 J= 1,9 

111 TS ( I I , J ) = TS ( I , J ) 

C 

CALL DC OS 2 ( C J, E J , S, KC J, KE J, KS ) S=DC 

CALL MUL.TA ( TL, S ,1 8, 18 , 1 8, KTL,KS ) 

IF (NAMES! .EG . 6h .OR. NAMEST .EQ. 6HN0STRS) GO TO 210 

CALL M'JLTA ( TS, S ,NRS T» 18 »1P» KTS ,KS ) 

210 CALL ATXBA1 ( S ,TL ,18, 18 ,KS,KTL ) 

RETURN 

C 

999 CALL ZZBOMB ( 5HSTE2 , NERROR) 

END 
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SUBROUTINE S.TE3 (C J, EJ , TMEM , TBEN, E ,ANU,NAMEK ,NAMEST ,S ,TL ,TS,NR ST * 

* KCJ ,KEJ,KS,KTL,KTS) 

DIMENSION CJ(KCJ,1) ,E J(KE J, 1 ) ,S ( KS ♦ 1 ) ,T L( KTL, 1 ) ,TS ( KTS , 1 ) 

DIMENSION CW(3,3), EW(3,3), Wl(18,lft), 

* IVine), 1V2( 16 ) » IV3U8), IV4<18) 

DATA KCW,KW1 / 3,18 / 

DATA IV1/ It 2 t 3 t 4, 5, 6, 7, 8, 9 1 1 Ot 1 1 ,1 2 1 13*1 4, 15 ♦ 16 , 17, 1 8/ , 

* 1V2/ I, 2 1 3 1 4, 5, 6,13,14,15,16,17,18,19,20,21,22,23,24/, 

* IV3/ 1, 2, 3, 4, 5, 6, 7, 8, 9 , 10, 1 1 , 1 2, 1 9,20,21 ,22,23 ,24/ , 

* 1 V 4/ 7, 6, 9,10, 11, 12,13,1 4 ,15, 16, 17, 18, 19, 20, 21, 22, 23 ,24/ 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT,.. 

C STIFFNESS MATRIX (SAME AS GLOBAL LOAD TRANSFORMATION MATRIX), 

C LOCAL LOAT TRANSFORMATION MATRIX (NOT YFT), 

C STRESS TRANSFORMATION MATPIX (NOT YFT), 

C FOR A COMBINED MEMBRANE-BENDING QUADRILATERAL PLATE ELEMENT WITH 
C UNRESTRAINED BOUNDARIES. 

C STIEFNESS MATRIX IS IN GLOBAL COORDINATE DIRECTIONS. 

C GLOBAL COORDINATE ORDER IS 

C (U,V,W,P,Q,R) JOINT 1, THEN JOINT 2, 3, 4. 

C WHERE U,V,W ARF TRANSLATIONS AND P,Q,P ARE POTATIONS. 

C GLOBAL LOAD TRANSFORMATION MATRIX PELATES LOADS AT QUAD VERTICES IN 
C GLOBAL COORDINATE DIRECTIONS TO DEFLECTIONS IN THE GLOBAL COORDINATE 
C DIRECTIONS. 

C ROW ORDFR IN GLOBAL LOAD TRANSFORMATION MATPIX IS 
C ( PU , RV, RW »MP , MR, MR ) JOINT 1, THEN JOINT 2,3,4. 

C WHERE P IS FORCE AND M IS MOMENT. 

C LOCAL LOAD TRANSFORMATION MATRIX RELATES LOADS AT QUAD VERTICES 
C IN LOCAL COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE 
C DIRECTION. 

C STRESS TRANSFORMATION MATRIX RELATES STRESS AT QUAD VERTICES IN LOCAL 
C COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE DIRECTION. 

C EULER ANGLF CONVENTION IS GLOBAL X,Y,Z PERMUTATION. 

C CALLS FORMA SUBROUTINES STF2 ,P EV ADD, ZZ BOMB . 

C DEVELOPED BY WA PENFIFLD, RL WDHLEN. FEBRUARY 1973. 

C LAST REVISION BY WA BENF1ELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS 

C CJ = INPUT MATPIX OF GLOBAL X,Y,Z COORDINATES AT QUAD JOINTS. 

C ROWS 1,2,3 CORRESPOND TO X,Y,Z COORDINATES. 

C COLS 1 ,2 ,3,4 CORRESPOND TO JOINTS 1,2, 3, 4. SIZE<3,4). 

C EJ = INPUT MATPIX OF EULER ANGLES (DFGREES) AT QUAD JOINTS. 

C ROWS 1,2,3 CORRESPOND TO GLOBAL X,Y,Z PERMUTATION. 

C COLS 1,2 ,3,4 CORRESPOND TO JOINTS 1, 2,3,4. S1ZE(3,4). 

C TMEM = INPUT EEEFC.T1VE MEMBRANE THICKNESS. 

C TBEN = INPUT EFFECTIVE BENDING THICKNESS. 

C E = INPUT YOUNGS MODULUS OF ELASTICITY. 

C ANU = INPUT POISSONS RATIO. (F/2G)-1. 

C NAMEK = INPUT TYPE OB STIF MATRIX WANT FD . 

C = Kit UStS A TRIANGLES, OVERLAP AVERAGE. 

C NAMEST = INPUT OPTION FOR S7RFSS TRANSFORMATION. 

C = 6H DR GHNOSTRS ,N0 STRESS TRANS CALCULATED. 

C S = OUTPUT STIFFNESS MATPIX (SAME AS GLOBAL LOAD TRANSFORMATION 

C MATRIX). SIZF(/4,24|. 

C TL = OUTPUT LOCAL LOAD TRANSFORMATION MATRIX. SIZE(24,24). 
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C TS = OUTPUT STRESS TRANSFORMATION MATRIX. SI 2 F (NPST , 24 ) . 

C NR ST - OUTPUT NUMPFR OF ROWS IN STRFSS TRANSFOPMAT ION MATRIX. 

C KCJ = INPUT ROW DIMENSION OF CJ IN CALLING PROGRAM. MIN=2. 

C KEJ = INPUT ROW D1MFNS ION OF EJ IN CALLING PROGRAM. MIN=3. 

C KS = INPUT ROW DIMENSION OF S IN CALLING PROGRAM. MIN=24. 

C KTL = INPUT ROW DIMENSION OF TL IN CALLING PROGRAM. MIN=24. 

C KTS = INPUT ROW DIMENSION DF TS IN CALLING PROGRAM. M1N=NRST. 

C 

C NERPCR EXPLANATION 

C I = SIZE LIMITATION EXCEEDED. 

C 2 = NAMEK IMPROPERLY DFF1NED. 

C 

NRST * 24 

NERROR=l 

IE (KS .LT. 24 .OR. KTL .LT. 2*, .OR. KTS .LT. NRST) GO TO 999 
DO 5 J=1 ,24 
DO 5 1=] ,24 
5 S(I*J) = P.0 

IF (NAMFK .EC. 6HM ) GO TO 110 

NERRCP=2 

GO TO 999 
C 

110 DO 200 1=1,3 

CWf 1,1 ) = CJ( 1,1 ) 

FW(I,1I = F J ( I , 1 ) 

CWf 1,2) = CJII,2) 

FW ( 1 , 2 ) = F J ( I , 2 ) 

CW ( 1,3) = C J( 1 ,3 ) 

200 EW ( 1 ,3 ) = F J ( 1 , 3 ) 

CALL STF2 ( CW,E /,TMEM ,T t FN, E ,ANU , NAMEK, NAMEST ,W 1 , TL , TS ,NRSTX, 

* KCW , KCW , KW 1 » KTL » KTS ) 

CALL PEVADD ( . ,W 1 , IV 1 , 1V1 , S, 18,18,24, ?4, 18,KS) 

DO 201 1=1,3 

cwd,i ) = c j ( : i ’ ) 

EW f 1 , 1 ) = EJf 1,1) 

C W ( 1 , 2 ) = C J ( I ,3 ' 

EW ( 1 , 2 ) = E J( I - 3 ) 

CW (1,3) = C J ( 1 , 4 ) 

201 EW 1 1 ,3 ) = EJ(I,4) 

CALL STF2 ( CW, l W, TM f M ,TDF N, F , AMU, NAMFK, NAMEST, W 1 ,TL,TS, NR STX, 

* KCW, KCW, KW1 , KTL, KTS ) 

CALL PEVAPO (•b,Wl,IV2,lV2,S, 18,18,24,24, 18, KS) 

DO 203 1=1,3 
CW( 1 ,1) = CJ( 1 ,1 ) 

E W ( I , 1 ) = F J ( I , 1 ) 

CW( I ,2 ) = C J ( I ,2 ) 

E W ( 1 , 2 ) = E J ( 1 , ? ) 

cw( i ,? ) = r j( l ,4) 

203 E W ( 1 ,3 ) = EJ(1,4) 

CALL STE2 ( CW , F W, TM FM ,T BE N, E , ANU , NAME K ,NAME ST ,W 1 ,TL , T S ,NRSTX , 

* KCW, KCW, KW 1, KTL, KTS ) 

CALL REVADD ( . 5, W 1 » IV3 » I V3 , S » 18,18,24,24, 18, KS) 

DC 205 1=1,3 
CW(1,1 ) = C J ( I , 2 ) 

e w ( i , i ) = ej(i t : ) 
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CW(I,2> = C J ( I *3 ) 

EW(1,2) = EJ<I,3) 

CVm,3) = CJ(I,4) 

205 FW(I,3I = EJ<1,4) 

CALL STF2 |CW,EW»TMEM»TBEN» E *ANU »NAMEK,NAMEST ,W1*TL,TS, NRSTX* 

* KCW,KCW,KW1,KTL,KTS) 

CALL REVADD < .5,W1 , IV4, IV4, S, 18,16,24,24, 18, KS) 

C 

DC 1 300 1,24 

DO 300 1=1,2*- 
TL ( I , J J = 0.0 
300 TS ( I , J ) = 0.0 
RETURN 
C 

999 CALL 2ZB0HB ( 4HSTF3 ,NERROR) 

END 
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SUBROUTINE STF3A (C J, EJ ,TH, G,NAMEK ,NAME ST,S ,TL ,TS ,NRST, 

* KCJ,KEJ,KS,KTL,KTS) 

DIMENSION CJ(KCJ,1), Ej(KEJ,l), S(KS,1), TL(KTL,1), TS(KTS,1) 

C 

C SUBROUTINE TO CALCULATE FINITE ELEMENT... 

C STIFFNESS MATRIX (SAMF AS GLOBAL LOAD TRANSFORMATION MATRIX)* 

C LOCAL LOAD TRANSFORMATION MATRIX, 

C STRESS TRANSFORMATION MATRIX, 

C FOP A PFCTANGULAP SHEAF PANEL ELEMENT WITH UNRESTRAINED BOUNDARIES. 

C STIFFNESS MATRIX IS IN GLOBAL COORDINATE DIRECTIONS. 

C GLOBAL COORDINATE ORDER IS 
C (U,V,W) JOINT 1, THEN JOINT 2, 3, 4. 

C WHERE U,V,W ARF TRANSLATIONS. 

C GLOBAL LOAD TRANSFORMATION MATRIX RELATES LOADS AT PANEL VERTICES IN 
C GLOBAL COORDINATE DIRECTIONS TO DEFLECTIONS IN THE GLOBAL COORDINATE 
C DIRECTIONS. 

C ROW OKDFR IN GLOBAL LOAD TRANSFORMATION MATRIX IS 
C (PU,PV,PW) JOINT 1, THEN JOINT 2, 3, 

C WHERE P IS FORCF. 

C LOCAL LOAD TRANSFORMATION MATRIX RELATES LOADS AT PANEL VERTICES IN 
C LOCAL COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE 
C DIRECTIONS. 

C ROW ORDFR IN LOCAL LOAD TRANSFORMATION MATRIX IS 
C FXl,PX2,PXr ,PXA, PY1 , PY2, PY3, PY4 

C WHERE P IS FORCE. X GOES FROM 1 TO 2 . Y GOES FROM 1 TO 4. 

C STRESS TRANSFORMATION MATRIX RELATES PANEL SHEAR STRESS (CONSTANT) IN 
C LOCAL COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORD DIRECTIONS. 
C EULER ANGLE CONVENTION IS GLOBAL X,Y,2 PERMUTATION. 

C CALLS FORMA SUBROUTINES ATXB A1 ,DC0S3C ,K3C1 ,MULTA ,ZZBCMB . 

C DEVELOPED BY RL WCHLEN. APRIL 1974. 

C LAST REVISION BY WA BENFIELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS 

C CJ = INPUT MATRIX OF GLOBAL X,Y,2 COORDINATES AT PANEL JOINTS. 

C ROWS 1,2,3 CORRESPOND TO X,Y,Z COORDINATES. 

C COLS 1,2, 3, 4 CORRESPOND TO JOINTS 1, 2,3,4. SIZE(3*4). 

C EJ = INPUT MATRIX OF EULER ANGLES (DEGPEES) AT PANEL JOINTS. 

C ROWS 1,2,3 CORRESPOND TO GLOBAL X,Y,Z PERMUTATION. 

C COLS 1,2 ,3,4 CORRESPOND TO JOINTS 1,2,3, 4. SIZEI3,4). 

C TH = INPUT PANEL THICKNESS. 

C G = INPUT SHEAR MODULUS OF ELASTICITY. 

C NAMEK = INPUT TYPE OF STIF MATRIX WANTED. 

C = M, USES K3C.1, 

C NAMEST = INPUT OPTION FOR STRESS TP AN SF OR NATION. 

C = 6H OF 6HNDSTR S ,NO STRESS TRANS CALCULATED. 

C S = OUTPUT STIFFNESS MATRIX (SAME AS GLOBAL LOAD TRANSFORMATION 

C MATRIX). S1ZF(12,12). 

C TL = OUTPUT LOCAL LOAD TRANSFORMATION MATRIX. SI ZE ( S , 1 2 ) , 

C TS = OUTPUT STRESS TRANSFORMATION MATRIX. S1ZE(1,12). 

C NPST = OUTPUT NUMBFP OF ROWS (1) IN STRESS TRANSFORMATION MATRIX. 

C KCJ =■ INPUT ROW DIMENSION OF Cj IN CALLING PROGRAM. 

C KEJ = INPUT »OW DIMENSION OF FJ IN CALLING PROGRAM. 

C KS = INPUT ROW DIMENSION OF S IN CALLING PROGRAM. MIN=12. 

C KTL = INPUT ROW DIMENSION OF TL IN CALLING PROGRAM. MIN=6 • 

C KTS = INPUT POW DIMENSION OF TS IN CALLING PROGRAM- M1N=1. 
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NEFROR EXPLANATION 

1 = SIZE LIMITATION EXCEEDED. 

2 = NAMFK IMPROPERLY DEFINED. 

NR ST = 1 

NERROR= 

IF IKS .LT. 1? .OR. KTL .LT. £ .OR. KT$ .LT. NR ST ) GO TO 999 
SL1 2 = SQPTUCJ(1*2)-CJ(1 ,1))**2 ♦ <C J( 2,2>-CJ(2» 1) 1**2 

* ♦ (CJ(3*2)— CJ(3»1)) **2 ) 

SL14 = SQ C TI ( C Ji 1 *4 )— C J (1*1)) **2 ■*- CC J ( 2,4)-C J C 2 t 1 ) )**2 

* ♦ (CJf 3,4)-CJI3,l)l**2> 

IF (NAMEK .EC. fcHKl ) GO TO 110 

N ERROR= 

GO TO 999 

110 CALL K3C1 (SL12»SL14,TH,G,TL,TS,KTL,KTS) 

C 

CALL DCOS ?C (CJ»EJ,S,KCJ,KEJ,KS) 

CALL MULTA ( TL , S ,E ,6 , 1 2, KTL »KS J 

IF (NAMEST .EC. frH .OR. NAMEST .EC. 6HN0STRS) GO TO 210 

CALL MULTA t TS, 5 ,NRST, fi , 1? ,KTS ,KS ) 

210 CALL ATXEA1 < S ,TL,E ,1 2,KS ,KTL ) 

RETURN 

C 

999 CALL ZZPOMF ( 6HSTF3A ,NERROR) 

END 


2 

TL=K 
S =DC 



TEGEOM 


SUBROUTINE TEGFOM <CJ,JM, VL.DV, KCJ, IFBAD ) 

DIMENSION C.*( KCJ * 1 ) * JM ( 1)* DV(1) 

DIMENSION R12(3) ,R13(3),R14(3) 

DATA EPS / l.E-5 / 

C 

C SUBROUTINE TO DETERMINE THE VOLUME AND VOLUMF CHANGE COEFFICIENTS OF 
C A TETRAHEDRON. 

C CALLS FORMA SUBROUTINES VCROSS ,VDOT . 

C DEVELOPED PY C S BCDLFY. FEBRUARY 19?4. 

C LAST REVISION PY R A PHIL1PPUS. AUGUST 1974. 

C SUBROUTINE ARGUMENTS 

C CJ = INPUT MATFIX OF JOINT CCOPDINATES. SIZE (3.8). 

C JM = INPUT VECTOR OF JOINTS DEFINING A TETRAHEDRON. SIZE (4). 

C VL = CUTPU1 VOLUME OF TETRAHEDRON DEFINED FY JM. 

C DV = OUTPUT VECTOR OF VOLUME CHANGE COEFFICIENTS. 

C KCJ = INPUT POW DIMENSION SIZE OF CJ IN CALLING PROGRAM. MIN = 3. 
C IFBAD = OUTPUT 

C = 0 * THE TETRAHEDRON VERTICIES ARE NOT NUMBERED ACCORDING 

C TO THE ESTABLISHED CONVENTION, OR LIE IN A PLANE. 

C 

J1 = JM ( 1 ) 

J? = JM(2 ) 

J3 = JM(3 > 

JA = JM (4 ) 

DO 5 1=1,2 

PI? (I) = CJ(ItJ?) - CJ(I.Jl) 

R13( I 1 = CJ( I , J3) - CJ(I.Jl) 

5 R14(I> = C J( I * J4 ) - CJ(I.JI) 

C 

CALL VCROSS < R 12 , P 1 3, DV < 10) , VAMAG , VPM AG ,VZMAG,SINAB ) 

CALL VDOT (DV(10),RI4,VCL,VAMAG,VBMAG,C0SAB) 

IF (VOL.LF.EPS) 1 FBAD=0 
VL = VOL/6. 

C 

CALL VC FOSS ( F 13 ,P 14, DV (4 ) , VAMAG,VPMAG ,VZMAG,SINAB) 

CALL VCROSS ( R 14,R12, DV (7 ) , VAMAG, VbMAG,VZMAG,SINAB ) 

DC 10 1 = 1 ,3 

10 OV(I) = -DV ( I +3 J - DV ( I +6 ) - DV < 1+9 ) 

DO IS 1=1 ,1? 

15 DVI1) = DVIII/6. 

C 

RETURN 

END 



TRNGL 


1/ 5 


SUBROUTINE TRNGL ( XYZ, JDCE ,EUL,NUTFL ,N J, 

* NU TMX, NUTKX, NUTBX, NUTLT, NUTST t 

★ W,T,$,KX,KJ,KE,KW) 

D1MFNSION <YZ(KX,1) ,JDOF(KJ,l ) ,EUL(KE,1 ),W(KW, I ),T(KW,1) ,S(KW,1 ) 
DIMENS I CM CJ(3,3>, EJ(3,3), IVH18) 

DATA NAMES /6HTRNGL /, NRW ,NRLT/18,18/, IBLNK/6H /, KCJ/3/ 

DATA «lT'.NCT/5*6/ 

C 

C SUBROUTINE TO CALCULATE CON OPTION) FINITE ELEMENT ... 

C MASS MA7R5LCFS AND IVECS (ON NUTMX), 

C STIFFNESS MATRICES (SAME AS GLOBAL LOAD TRANSFORMATION MATRICES) 

C AND IVECS (ON NUTKX ) * 

C UNIT LOAD BUCKLING MATRICES AND IVECS (ON NUTBX), (NOT YET) 

C LOCAL LOV5 TRANSFORMATION MATRICES AND IVECS (ON NUTLT), 

C STRE SS r P sNSFORMATION MATRICES AND IVECS (ON NUTST), 

C FOR COREINEC MEMBRANE-PENDING IRIANGLF PLATE ELEMENTS. 

C MASS, STIFFNESS, BUCKLING MATRICES ARE IN GLOBAL COORDINATE 
C DIRECTIONS * 

C GLOBAL COORDINATE ORDER IS 

C (U,V,V(,P,!V,R) JOINT I , THEN JOINT 2, 3. 

C WHERE U,V,W APE TRANSLATIONS AND P,0,P APE ROTATIONS. 

C IVEC GIVES ELEMENT DO r INTO GLOBAL DOF. EXAMPLES... 

C IVFC(6i-.03~ PLACES ELEMENT DCF 6 INTO GLOBAL DCF 834. 

C 1 VFC ( 3 : =0 OMITS ELEMENT DOF 3 FROM GLOBAL DOF. THIS CONSTRAINS 

C ELEMENT DCF 3 TO ZERO MOTION. 

C GLOBAL LOAD TRANSFORMATION MATRIX RELATES LOADS AT TRNGL VERTICES IN 
C GLOPAL CDCRDINATF DIRECTIONS TO DEFLECTIONS IN THE GLOBAL COORDINATE 
C DIRECTIONS. 

C ROW ORDER IN GLOBAL LOAD TRANSFORMATION MATRIX IS 
C ( PU, PV , PW ,MP»MO?MR ' JOINT 1, THEN JOINT 2,3. 

C WHERE P IS FORGE AND M IS MOMENT. 

C LOCAL LOAD TRANSFORMATION MATRIX RELATFS LOADS AT TRNGL VERTICES IN 
C LOCAL COORDINATE SYSTEM TO DEFLECTIONS IN THE GLOBAL COORDINATE 
C DIRECTIONS. 

C ROW ORDER IN LOCAL LOAD TRANSFORMATION MATRIX IS 
C (PX,RY,MZ) JOINT 1 THEN 2,3, NEXT 

C (PZ,MX,MY> JOINT I THEN 2,3. 

C WHERF P IS FORCE AMD V IS MOMENT. 

C STRESS TRANSFORMATION .-A TP. IX RELATES STRESS AT TRNGL VERTICES IN LOCAL 
r. COORDINATE SYSTEM TO' DEFLECTIONS IN THE GLOBAL COORDINATE DIRECTIONS. 

C ROW CRDFP. IN STRE G TRANSFORMATION MATRIX IS 
C ( SIGMA— X , SIGMA,— Y, TAU— XY) FOR (Z=TBEN/2) AT JOINT 1, 

C THEN JOINT 2,3. 

C ( SIGMA-X, c IGMA-Y,TAU-XY) FOR (Z=-TEEN/2) AT JOINT 1, 

C THEN JOINT 2,3. 

C WHERE SIC . IS NORMAL STRESS AND TAU IS SHEAR STRFSS. 

C DATA AR *ANGEMEf!T ON NUTMX, NUTKX, NUTBX, NUTLT, NUTST FOR EACH 
C FINITE ELEMENT IS ( W-M ,K ,b »LT * ST ) 

C WRITE (NUTMX) N».ME W ,N EL ,NP , NC , NA.MFL , ( IB LNK, 1 = 1,5) , 

C ((W(1,J),1=1»NR),J=1,NC)»( IVEC ( I ) »I-1*NC ) 

C CALLS FORMA SUBROUTINES M AS 2 , P AGEHD , STF2, ZZBOME. 

C DEVELOPED c ^ WA BFMEIELD, CS BODLEY, Rl WOHLEN. FEBRUARY 1973. 

C LAST PF7TMPN FY RL WOHLEN. MAY 1976. 

C 

C A* V V* fr**:?**** $:$:**■*:#:*##***#♦*#*#** 
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INPUT DATA READ IN THIS SUBROUTINE FROM NUTEL. 
READ FROM CARDS. 

NAMFM,NAMEK,NAMELT,NAMFST,NAME& 

RO» E » ANU 

TMASC,TMEMC,TBENC 
20 NEL»J1,J2 T J3 *TM/ ^ VtTMEM V» TBEN V 
IF (J1 . EQ. 0) RETURN 
GO TO 20 


IF NUTEL = 5, DATA IS 

FORMAT (5(A6,<tX) 
FORMAT (3(5X,E10)) 
FORMAT (3(BX,E10I) 
FORMAT (AI5,3E10) 


DEFINITION OF INPUT VARIABLES. 

NAMEM = TYPE OF MASS MATRIX WANTED. 
= Ml, DIAGONAL LUMPED. 

= M2, CONSISTENT. 


NAMEK 


NAMELT 

NAME ST 

NAMEP 

RO 

E 

ANU 

TMASC 

TMASV 

TMEMC 

TMEMV 

TBENC 

TBENV 

NEL 

J1 

J2 

J3 


- 6H OR GHNCMaSS, NO MASS MATRIX CALCULATED. 

TYPE OF STIFFNESS MATRIX WANTED. 

= Kl, QUADRATIC DISPLACEMENT FOR MEMBRANE, CUBIC 
DISPLACEMENT FOR BENDING. 

= 6H OR 6HN0STIF, NO STIFFNESS MATRIX CALCULATED. 

IDENTIFICATION NAME FOR LOAD TRANSFORMATION MATRICES. 

= 6 H OR 6HN0LCAD, NC LOAD TRANSFORMATIONS CALCULATED. 

IDENTIFICATION NAME FOP STRESS TRANSFORMATION MATRICES. 

= 6H OR 6HN0STP S » NO STRFSS TRANSFORMATIONS CALCULATED. 

TYPE OF BUCKLING MATRIX WANTED. 

= 6H OR 6HN0BUCK, NO BUCKLING MATRIX CALCULATED. 


MASS DENSITY. 
YOUNGS MODULUS 
POISSONS RATIO. 
EFFECTIVE MASS 
EFFECTIVE MASS 
1 F .LF. 0., TMA 
EFFECTIVE MEM PR 
FFFECTIVE MEM BP 
IF .LF. 0 . , TME 
EFFECTIVE EENDI 
EFFECTIVE EENDI 
IF .LE. 0., TBE 
FINITE ELEMENT 
CALCULATIONS. W 
JOINT NUMBER AT 
JOINT NUMFEP AT 
JOINT NUMBER AT 


OF ELASTICITY. 

( E/2G )— I . 

THICKNESS, 

ThICKNESS, 

SC IS USED. 

ANF THICKNESS, 

ANE THICKNESS, 

MC IS USED. 

NG THICKNESS, 

NG THICKNESS, 

NC IS USED. 

NUMEER. FOR REFERENCE ONLY, 
RITTEN ON NUTMX, ETC. 
TRIANGLE VERTEX 1. 

TRIANGLE VERTEX 2. 

TRIANGLE VERTEX 3. 


I CONSTANT ) . 
(VARIABLE). 

(CONSTANT). 

(VARIABLE). 

(CONSTANT). 
(VARIABLE ). 


NOT USED IN 


EXPLANATION OF INPUT FORMATS. NUMBER INDICATES CARD COLUMNS USED. 

I - 1NTFGFR DATA, RIGHT ADJUSTED. 

E = DECIMAL POINT DATA, ANYWHERE IN FIELD. EXPONENT RIGHT ADJUSTED. 

X = CARD COLUMNS SKIPPED. 

#*** #3}: if:#*:*## 3 ^* + ««*** ********** ***$***************** * 


SUPPOUTINf A°CUMFNT$ (ALL INPUT) 

XY2 - MATRIX OF JOINT GLOBAL X,Y,Z LOCATIONS. ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 
X,Y,Z LOCATIONS RESPECTIVELY. SIZE(MJ,3). 

JDOF - MATRIX OF JOINT GLOBAL DEGREES OF FREEDOM. ROWS CORRESPOND 
TP JOINT NUMBERS. COLUMNS 1,2,3 CORRESPOND TO THE JOINT 
TRANSLATION DOFS AND COLUMNS CORRESPOND TO THE JOINT 
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C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


EUL 


NUTEL 

NJ 

NUTMX 


NUTKX 


NUTBX 


NUTLT 


NUTST 


W 

T 

S 

KX 

KJ 

KE 

KW 


ROTATION DOES. SIZE(NJ,6). 

MATRIX OF JOINT EULER ANGLES (DEGREES). ROWS CORRESPOND 
TO JOINT NUMBERS. COLUMNS 1*2,3 CORRESPOND TO THE 
GLOBAL X,Y,Z PERMUTATION. SIZE(NJ,3). 

LOGICAL NUMBER OF TAPE CONTAINING ELEMENT INPUT DATA FOR 
THIS SUBROUTINE. IF NUTEL = 5, DATA IS READ FROM CARDS. 
NUMBER OF JOINTS OR ROWS IN MATRICES (XYZ), (JDOF), (EUL). 
LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 
MASS MATP1CES AND IVECS ARE OUTPUT. 

NUTMX MAY BE ZERO IF MASS MATRIX IS NOT FORMED. 

USES FORTRAN READ, WRITE. 

LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 
STIFFNESS MATRICES (SAME AS GLOBAL LOADS TRANSFORMATION 
MATRICES) AND IVECS ARE OUTPUT. 

NUTKX MAY BE ZERO IF STIFFNESS MATRIX IS NOT FORMED. 

USES FORTRAN READ, WRITF. 

LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT UNIT LOAD 
BUCKLING MATRICES AND IVECS ARE OUTPUT. 

NUTBX MAY BE ZERO IF BUCKLING MATRIX IS NOT FORMED. 

USES FORTRAN READ, WFITE. 

LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT LOCAL 
LOAD TRANSFORMATION MATRICES AND IVECS ARE OUTPUT. 

NUTLT MAY BE ZERO IF LOAD TRANSFORMATIONS ARE NOT FOkMED. 
USES FORTRAN READ, WRITE. 

LOGICAL NUMBER OF UTILITY TAPF ON WHICH ELEMENT 
STRESS TF AN S FOR MAT I ON MA1RICES AND IVECS ARE OUTPUT. 

NUT ST MAY BE ZERO IF STRESS TRANSFORMATIONS ARE NOT FORMED. 
USES FORTRAN READ, WRITE. 

MATRIX WORK SPACE. MIN SIZE(18,18). 

MATRIX WORK SPACE. MIN S1ZF(1B,18). 

MATRIX WOPK SPACE. MIN SIZE ( 1 B , 18 ) . 

POW DIMENSION OF XYZ IN CALLING PROGRAM. 

ROW DIMENSION OF JDOF IN CALLING PROGRAM. 

ROW DIMENSION OF FUL IN CALLING PROGRAM. 

ROW DIMENSION OF W, T, AND S IN CALLING PROGRAM. MIN=18. 


NERPOR EXPLANATION 

1 = JOINT NUMBER GREATER THAN NUMBER OF JOINTS. 

2 = MASS MATRIX FORMED, NUTMX .LE. ZERO. 

3 = ST1FFNFSS MATRIX FORMFD, NUTKX .LE. ZERO. 

4 = LT MATRIX FORMED, NUTLT .LE. ZERO. 

5 = ST MATRIX FORMED, NUTST .LE. ZERO. 


1001 FORMAT 

1002 FORMAT 

1003 FORMAT 

2001 FORMAT 
* 

2002 FORMAT 
* 

2003 FORMAT 
* 

* 

* 

* 


(5<A6,*X) ) 

(3(BX,E10.P) ) 

(41B,3F10.u) 

<//??X ^9H]NPUT LATA FOR COMBINED MEMBRANE-BENDING TRIANGLE 
1EH PLATF FLFMENTS ) 

(//26X A9HINPUT DATA FOR COMBINED MEMBRANE -BEN DING TRIANGLE 
27H PLATE ELEMENTS (CONTINUED)) 

(/ 13X7HMASS = A6, 1 3X7HST I F = A6, 6X1 3HL0AD TRANS = A6, 

3X1 EE'S T PE SS TRANS = A6, 3X1 1HPUCKLING = A6, 

/ TbXAHPP = E10.3, 13X3HE * T10.3, 

/ 10X C, HT (MASS ) = E10.3, 12X*PNU - E10.3, 

/ 3?X1 3HT (MEMBRANE ) ~ E10.3, 



noon on 
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* / 33X1 2HT (BENDING ) = E10.3, 

* //18X 7HELEMENT 5X 7HJ0INT 1 5X 7HJ01NT 2 5X 7H JOINT 3 

* 5X 7HT(M AS S ) 6X 11HT (MEMBRAN E ) 5X 10HT (BENDING ) 

* /18X 6HNUMBEP 36X 3(5X 10H ( VAR IABLE ) ) ) 

2004 FORMAT (IPX 4(15»7X)»3(F10.3»5X) ) 

2005 FORMAT (IPX 4(15, 7X) ) 

READ AND WRITE FINITE ELEMENT DATA. 

NLINE = c 
CALL PAGE HO 
WRITF (NCi,?001) 

RFAD (NIJT C « ,1001) NAM EM,N AMEK ,NAMELT»NAME$T »NAME6 
READ (N0TFLvl002) RO,t,ANU 
RFAD (NLiTELt 1002 I TMASC »TMEMC , TBENC 
WRITE (NOT, 2003) NAME M,NAMEK, NAMELT,NAMEST,NAMEB, 

* RO,E ,TMASC ,ANU,TMEMC, TBENC 
20 READ ( NUT EL , 1 003 ) NEL , J 1 , J2 , J3 , TM AS V,TM EMV,TB ENV 

NO THIK = 1 

IF (TMaSV.LF.O. .AND. TKEMV.LE.O. .AND. TBENV.LE.O.) NO TH1K=0 
IF (J1 .LE. 0) RETURN 
NLINE = NLINE + 1 
IF (NLINE .LE. 42) GO TO 30 
CALL PAGEKO 
WRITE (NOT, 2002) 

WRITE (NOT, 2003) NAME M, NAME K, NAME LT , NAM EFT *NAME6* 

* RC,E , TMASC, ANU,TMEMC,TFENC 
NLINE = 0 

30 IF (NO THIK. EC. 1 ) 

♦WRITE ( NOT , 2004 ) NF L, J1 , J 2, J3 ,TMASV ,TME MV.TPENV 
IF (NO THIK. EC. 0) WRITE (NOT, 2005) NEL,J1*J2,J3 

NERR0R=1 

IE (JJ .GT. NJ .OR. J 2 .GT. NJ .OR. J3 .GT. NJ) GO TO 999 

SET THICKNESSES.' 

TMAS = TMASC 
TMEM = TM EMC 
TEEN = TBENC 

IF (TMASV.GT.,0.) TMAS=TMASV 

;f (tmfmv.gt.o. ) tmfm=tmemv 

,r (TPENV.GT.O. ) TEEN=TPENV 

FORM FINITE ELEMENT COORDINATE LOCATIONS, EULER ANGLES, REVADD IVEC. 
DO 42 I =1 ,3 
CJ(I,I) = XYZ ( J1 ♦ I ) 

C J ( I » 2 ) = XY2 ( J2 , I ) 

C J ( 1 ,3 ) = X.Y2 ( J3 , I ) 

E J ( 1 , 1 ) - FUL ( J1 ,1) 

E J ( 1 , ? J = EUL ( J2 » I ) 

42 E J ( I , 3 ) - EUL ( J3 , I ) 

DU 4^ 1=1 ,6 
I VI (I) = JOPF ( J1 , 1 ) 

IV1 (I+fc ) = JDOFl J2,l ) 

44 1V1U + 12) =■ JDOF (J3,I ) 

C 

C FORM MASS MATRIX (W). 
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IF (NAMFM .EQ. 6H .OR. NAMEM .EO. 6HN0MASS) GO TO 110 

CALL MAS2 ( C J, E J , TM AS ,P 0, NAMEM ,W , T, S, KCJ, KCJ , KW ,KW, KW ) 

NERR0R=2 

IF t NUT MX .LE. 0) GO TO 999 

WRITE (NUTMX) MAMFM t N EL ,NRW ,NRW,NAMFL , ( IBLNK, 1 = 1 , 5 ) , 

* ((W(1,J ),I=1,NRW),J~1 ,NPW),(lVltl>,l=l,NRW) 

C 

C FORM STIFFNESS MATRIX (W ), LOCAL LOAD TRANSFORMATION MATRIX (T), 

C STRESS TRANSFORMATION MATRIX (S>- 

110 IF (NAMEK .FQ. 6H .OR. NAMFK .FO. 6HN0STIF) GO TO 20 

CALL STF2 ( C J , E J , TM EM ,T6EN, F , ANU t NAME K,NAMEST,W ,T, S ,NRST t 

* KCJ,kCJ,KW ,KW,KW ) 

NERRCR=3 

IF (NUTKX .LE. 01 GO TO 999 

WPITF (NUTKX) NAMEK ,NEL ,NPK ,NRW,NAMEL ,( IPLNK, 1=1 , 5 ) , 

* UW(I,J),I=1,NRW),J=1,NRW),IIV1 (I),I=1,NRW) 

IF (NAMFLT .EC. fcH .OR. NAMFLT .EC. 6HN0L0AD ) GO TO 115 

NERR0R=9 

IF (NUTLT .LF. C) GO TO 999 

WRITE (NUTLT) NAM ELT, NE L , NR LT ,NRW,N AMEL * ( IP LNK , 1= 1 , 5 ) * 

* ( (T (I ,J ),I=1,NRLT) ,J=1,NPW) »( IVl(I) ,1=1 ,NRW) 

115 IF (NAMFST .FC. 6H .OR. NAMEST .tt. 6HN0STRS) GO TO 20 

NERR0R=5 

IF ( NU T c T .Lr. 0) GO TO 999 

WRITE (MUTST) NAME ST, NEL, NRST ,NRW,N AMEL , ( IBLNK , 1= 1 , 5 ) * 

* ( (SCI » J ), 1=1, NR ST) * J=1*NRW ),(IV1(1) , 1 =1 ,NRW ) 

GO TO 2C 

C 

999 CALL ZZBOMB (6HTRNGL ,NEkROR) 

END 



